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Abstract. Using a non-perturbative metliod developed in a previous article (aaper 11) we investigate the tails 
of the probability distribution V{pr) of the overdensity within spherical cells. Since our approach is based on a 
steepest-descent approximation which should yield exact results in the limit of rare events it applies to all values 
of the rms linear density fluctuation a, from the quasi-linear up to the highly non-linear regime. First, we derive 
the low-density tail of the pdf. We show that it agrees with perturbative results when the latter are finite up 
to the first subleading term, that is for linear power-spectra P{k) oc fc" with —3 < n < —1. Over the range 
— 1 < n < 1 some shell-crossing occurs (which leads to the break-up of perturbative approaches) but this does not 
invalidate our approach. In particular, we explain that we can still obtain an approximation for the low-density 
tail of the pdf. This feature also clearly shows that perturbative results should be viewed with caution even when 
they are finite. We point out that our results can be recovered by a simple spherical model (this is related to the 
spherical symmetry of our problem). On the other hand, we show that this low-density tail cannot be derived from 
the stable-clustering ansatz in the regime a S> 1 since it involves underdense regions which are still expanding. 
Second, turning to high-density regions we explain that a naive study of the radial spherical dynamics fails. 
Indeed, a violent radial-orbit instability leads to a fast relaxation of collapsed halos (over one dynamical time) 
towards a roughly isotropic equilibrium velocity distribution. Then, the transverse velocity dispersion stabilizes 
the density profile so that almost spherical halos obey the stable-clustering ansatz for —3 < n < 1. We again 
find that our results for the high-density tail of the pdf agree with a simple spherical model (which takes into 
account virialization) . Moreover, they are consistent with the stable-clustering ansatz in the non-linear regime. 
Besides, our approach justifies the large-mass cutoff of the Press-Schechter mass function (although the various 
normalization parameters should be modified). Finally, we note that for cr ^ 1 an intermediate region of moderate 
density fiuctuations appears which calls for new non-perturbative tools. 
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1. Introduction 



In usual cosmological scenarios large-scale structures in 
the universe have formed through the growth of small ini- 
tial density fluctuations by gravitational instability (e.g., 
Peebles (1980)). Moreover, in most cases of cosmological 
interest the amplitude of these density fluctuations in- 
creases at smaller scales, as in the standard CDM model. 
This leads to a hierarchical scenario of structure forma- 
tion where smaller scales enter the non- linear regime first, 
building small virialized objects which later become part 
of increasingly large and massive structures. These halos 
will produce galaxies or clusters of galaxies (depending on 
non-gravitational processes like colhsional cooling) which 
build a complex network among large voids. Thus, in or- 
der to describe the non-linear structures we observe in the 
present universe it is of great interest to understand the 
non-linear evolution of the density field. 



Unfortunately, this is a rather difficult task since this 
non-linear regime (i.e. small scales or rare large density 
fluctuations on large scales) cannot be described by per- 
turbative methods (except large voids for n < —1). In 
fact, very few results have been obtained so far in this do- 
main since most rigorous approaches to the dynamics of 
gravitational clustering have relied on perturbative meth- 
ods. Therefore, they were restricted to the early stages of 
the non-linear evolution. In this paper, we make use of 
a non-per turbative method developed in a previous paper 
( paper 11 ) to tackle this fully non-linear regime. More pre- 
cisely, since this approach is based on a steepest-descent 
approximation we investigate the regime of rare events 
where it should yield asymptotically exact results. Thus, 
we consider the very high density and low density tails of 
the probability distribution function (pdf) V{pr) of the 
overdensity within spherical cells. Then, this approach ap- 
plies to any value of the rms linear density fluctuation a 
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SO that it describes all regimes, from linear to highly non- 
linear scales. However, while in the quasi-linear regime 
it is able to predict the whole pdf (this was investigated 
in a specific manner in paper II) in the highly non- linear 
regime it only provides the rare-event tails of the pdf. 

This article is organized as follows. First, in Sect.p] we 
recall the path-integral formalism developed in [paper 11 
which provides an explicit expression for the pdiV{pR) in 
terms of the initial conditions. Then, in Sect.H we inves- 
tigate the very low density tail of the pdf. We first derive 
the s add le-point of the actio n wh ich governs this regime in 
Sect.SJ and we give in Sect. |3.5| the lo w d ensity tail of the 
pdf this implies. We compare in Sect JS.q our results with 
perturbative methods, a simple spherical model and the 
usual stable-clustering ansatz. Next, in Sect.Q we turn to 
the high- den sity tail of the pdf. We first derive in Sect. 4.1 
and Sect. 4.2 a naive spherical saddle-point but we show in 
Sect. 4.3 that virialization processes are sufficiently violent 
to significantly modify this simple approach. Then, we give 
in Sect. 4.5 the high-density tail of the pdf obtained by a 
more careful study. We finally c omp are this result with the 
stable-clustering ansatz in Sect .[4.7| and with the standard 
Press-Schechter mass function (Press & Schechter (1974)) 



in Sect. 4.8 



2. Action S[5l] 

In this section, we introduce our n otations a nd we briefiy 
recall the formalism developed in paper II which allows 
us to evaluate the pdf V{5r). More precisely, we investi- 
gate the statistical properties of the overdensity pR within 
spherical cells of comoving radius i?, volume V: 

Pr = 1 + Sr with Sr = / —— (5(x). (1) 

Jv ^ 

Here (5(x) is the non-linear density contrast at the comov- 
ing coordinate x, at the time of interest. We investigate in 
this article the case of Gaussian initial conditions, so that 
the statistical properties of the random linear density field 
^l(x) are fully defined by the two-point correlation: 



Ai(xi,X2)=(5L(Xi)(5i(x2)). 



(2) 



The kernel Al is symmetric, homogeneous and isotropic: 
Al(xi,X2) = Aidxi— X2I). Moreover, the matrix A^ can 
be inverted and its inverse A^ is positive definite since 
we have: 



Sr.AV\Sr = / dk 



I^L(k)|^ 
P{k) 



for real fields (5l(x), see paper II. Here we introduced the 
short-hand notation: 



h-A-'.h 



/ dxidx2 /i(xi).A^^(xi,X2)./2(x2) 



(4) 



for any real fields /i and /2. Moreover, the power- 
spectrum P{k) of the linear density contrast used in eq.(||) 
is given by: 



((5L(ki)<5L(k2)) = P(fci)^Z3(ki 



where we defined the Fourier transform of the linear den- 
sity field by: 



5l(x) 



dke 



ik.x 



^L(k). 



(6) 



Finally, it is convenient to introduce the usual rms linear 
density fluctuation a{R) in a cell of radius R: 



a'{R) 



{SIr) 



dxi dx2 



(7) 



where Sl.r. is the mean linear density contrast over the 
volume V as in eq.(|]). 

The previous expressions describe the statistical prop- 
erties of the initial conditions, through the linearly- 
evolved density field (5l(x) (see also ^aper I ). Then, if 
we could explicitly solve the coUisionless Boltzmann equa- 
tion which governs the gravitational dynamics (coupled to 
the Poisson equation) this would provide the statistics of 
the actual non- linear density field 5(x). In any case, we 
can always write a closed formal expression for the pdf 
'P{Pr) of the overdensity pR, following the method pre- 
sented in paper II, To do so, it is convenient to introduce 



the Laplace transform "0(2/) of the pdf, given by: 



V-Cy) 



== lf>-VPR\ = 



dpR e-yP-- VipR) 



(8) 



since pR > 0. Here, the symbol (..) expresses the average 
over the initial conditions. Then, the pdf can be recovered 
from 4>{y) through the inverse Laplace transform: 



V{pr) = 



-\-ioo 



2m 



(9) 



As described in paper 11 , we can write an explicit expres- 
sion for the average over the initial conditions which ap- 
pears in eq.(|^). Indeed, since these initial conditions can 
be fully described by the hnear growing mode 5l(x) (see 
paper 1) which is a Gaussian random field we can write: 



V'(y) - (DetA, 



,1/2 



:d^L(x)] 



p-ypj?.[i5t]-^<5i--Aj^\(5£, 



(10) 



where the inverse kernel A^ was introduced in eq.(0). 
Here the functional pr[5l\ is the non- linear overdensity 
PR produced by the linear density field 5^. 

As in paper H, it is actually convenient to introduce 
the rescaled generating function ipiv) defined by: 



(3) i,{y) = i,{va^{R)) 



(11) 



Indeed, this allows us to factorize the amplitude of the 
power-spectrum P{k) in the "action" S[6l\ which charac- 
terizes our system. Thus, we can now write eq.(|lO|) as: 



^{y) ^ {miAl'Y'^ 



[d5L(x)] e 



-S[5z,]/<T^(fl,) 



where we introduced the action 5[(5l]: 
a\R) 



(5) S[5l\ = y pr[6l\ 



Sl.Ai\Sl 



(12) 



(13) 
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The action S[6l] is independent of the normalization of 
the linear power-spectrum P{k) since A^ oc cr^, see eg. (ji\) . 
In paper II the change of variable y -^ y/o"^ in eg. (pi]) 
was crucial since it allowed us to show that the steepest- 
descent method was asymptotically exact in the limit a -^ 
0. Indeed, it is clear that in this guasi-linear limit the path- 
integral in eg. (O) is given by the global minimum of the 
action S. However, in this article we no longer study this 
guasi-linear limit. Indeed, we investigate the regime of rare 
events, i.e. the limits pr ^ (rare "voids") or pn -^ -|-cx) 
(very high overdensities) , and we take a to be finite. Thus, 
our study applies both to the linear and non-linear regimes 
since the value of a is actually irrelevant. Then, the change 
of variable introduced in eg. ( |ll|) is no longer essential and 
we could directly work with ^p{y). Nevertheless, it is still 
convenient to work with the rescaled generating function 
ij^iy) since it allows us to get rid of the amplitude of the 
rms density fluctuations which plays no key role in the 
physics we investigate here. Besides, it allows us to use 
the same expressions as in paper 11 . 



3. Rare underdensities 

We first consider the statistical properties of very rare 
underdensities (i.e. "voids"), using the tools laid out in 
the previous section. Indeed, this is much easier than the 
study of large overdensities which involves shell-crossing 
in an essential way. Hence it is convenient to start with 
voids to recall the basic ideas behind the steepest-descent 
method introduced in jpaper II since it only reguires minor 
modifications to be applied to underdensities. We shall 
investigate the high-density tail of the pdiVipn) in a next 
section. 

In the following we restrict ourselves to the case of a 
critical-density universe. Then, the spherical solution of 
the coUisionless Boltzmann eguation is explicitly known, 
as long as there is no shell-crossing. 



3.1. Spherical saddle-point 

Thus, our goal is now to evaluate the path-integral ( [l2| ) in 
the limit pfj —>■ in order to derive the low-density tail of 
the pdiV{pB)- To this order, we can try a steepest-descent 
approximation, in the spirit of the calculation performed 
Since pr[5l] 



in paper H 



1 + ^fl[<5L] the action S [6l] 



is directly related to the action studied in paper II and 



we can use the results obtained in that previous work. In 
particular, if there is no shell-crossing we know that the 
action S[Sl] admits a spherically symmetric saddle-point 
(5l(x) given by: 



U,_Ri 



K[SL,R,]<yHRL)/<y^{R) 



1 - ^', [Sl,r,] 



L,Rl\ j^5l,r 



i^ <j(Rl) dR 



(Rl) 



together with: 



(5i,(x) =6l,r^ / -Tj- 



Vl 



dx^ Al(x,xO 
Vl cj^{Rl) 



(14) 



(15) 



Here the Lagrangian comoving radius Rl is such that the 
matter enclosed within this volume Vl in the primordial 
universe (i.e. Ml — 47r/3 J^Rl) ends up within the ra- 
dius R in the actual non-linear density field. Moreover, 
for such spherically symmetric initial conditions the ac- 
tual non-linear overdensity pR only depends on the linear 
density contrast 5l,Ri^ within the radius Rl through the 
function J- p. This function is the usual spherical solution 
of the dynamics. It can be expressed in terms of hyperbolic 



functions as (for r2,ji = 1, see Peebles (1980)) 



^p{5l) = 7T 



9 {s\nhrj — r]Y 



2 (cosh 77 -1)3 



'^^ ^ ~20 [^(^i'^l^'?^'^)]^'''' 



(16) 



Note that since we study pR — 1 + 5r rather than 5r 
the function Tp de fined in e g.([l6[) differs from the usual 
function T (used in paper II ) by a factor -1-1. Besides, the 
non-linear guantities pR and R are related to the linear 
variables 5l,Ri, and Rl by: 



PR = ^p [5l.,Rl\ 



Rl = PR R 



(17) 



Then, the implicit eg.(|l4|) defines the normalization Sl,r,^ 
of the spherical saddle-point associated to a given value of 
y, while eg.(|l5[) provides the density profile of this saddle- 
point. The cumulative linear density profile Sl,r' /^l.Rl of 



this spherical saddle-point is displayed in Fig. 1 in paper II 



We shall simply recall here that for a power-law linear 
power-spectrum P{k) ex fc" we have the asymptotic be- 
haviours: 



R'r 



and 



R'r 



ol,r'^ 

SlMl 



^„ (l-n)(3-n) 



JL.R' 



^L,Rl 



2" 



(l-n)(3-n) (Ri^ 

R'r 



(18) 



(19) 



We display in Fig.^ the cumulative non-linear overden- 
sity profile pri. It is obtained from the relation pRi = 
Tp[Sl,r' ] which apphes to all radn {R',R'^), where i?^ 
is the linear Lagrangian radius associated with the actual 
non- linear radius R' (assuming there is no shell-crossing). 
Of course, as for the linear density contrast Sl^r/ we 
find that pR' remains finite for small R', of the order oi pr, 
while it goes to unity (i.e. the density contrast vanishes) 
at large radius R'. Besides, since we now probe the highly 
non-linear regime we must pay attention to shell-crossing. 
Here, we must recall that the previous results were derived 
with the assumption that there is no shell-crossing, so that 
the simple mapping (|l^) is valid. However, it is clear in 
Fig.|l|that the profile obtained for the case n — 0.5 leads to 
shell-crossing at R' ~ R. The condition which expresses 
that shell-crossing appears is d_R'/di?^ = 0. Using this 
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Fig. 1. Cumulative non-linear overdensity profile pm of 
the spherical saddle-point. The solid line corresponds to 
n = —1.5, the dashed-line to n = —0.5 and the dot-dashed 
line to n = 0.5. For the curve n = 0.5 we clearly see that 
shell-crossing occurs at B! ~ R. 



constraint in the system ( |17| ) written for an arbitrary pair 
(R' ,R'^) we obtain: 



shell-crossing 



dlnJF„ 



dln(-<5i,fl;) 



>3 



dlni?; 



d\n{-5LM,' 



(20) 



Thus, we need the behaviour of the function TpiS^) for 
large negative S^- From eq.(p6|) we obtain at large 77: 



^p^Sl) 



20 a 



-3/2 



(21) 



Substituting this result into eq.(EO[) we get the condition: 



shell-crossing 



dln(-,5i^ 



R'r 



din R'r 



< -2. 



(22) 



Then, from the asymptotic behaviour (|l9|) we find that 
for n > —1 some shell-crossing occurs at scales R' ^ R for 
very underdense saddle-points. In Fig.|lJ the curve shown 
for the case n = —0.5 does not exhibit any shell-crossing 
yet because the normalization 5^^^ is not large enough 
but we can check numerically that for a more negative 
value of 5l,Ri, [pr ^ 10-^) some shell-crossing indeed 
appears. For n > some shell-crossing also occurs at R' ^ 
R since for such power-spectra the local density contrast 
^l(x) changes sign at the radius Rl- Then the derivative 
of the cumulative density contrast Sj^ j^/ is discontinuous 
at Rl and dlii{-SL.,R'J/dliiR'^ < -3 for R'^ > Rl. 

These results imply that the spherical saddle-point ob- 
tained in eq.(n) and eq.dlS) is not correct for large neg- 
ative Sl,Rl if '^ > ~1- However, since shell-crossing only 
appears over a limited range of radii of the order of Rl 
we can expect the previous state (5l(x) to be a reasonable 
approximation to the exact spherically symmetric saddle- 
point. Moreover, for — 1 < n < shell-crossing only occurs 
for R'j^ > Rl so that the radial profile (|l^) is correct for 
R'l < Rl and it should only be modified at i?^ > Rl- 



In order to derive the exact saddle-point we need to ex- 
plicitly take into account shell-crossing which makes the 
calculation much more difficult, since we do not know the 
exact functional /^^[(^l] in this regime. However, because of 
the spherical symmetry of the problem which we consider 
here, the saddle-point should remain spherically symmet- 
ric. In the following, we shall use the linear state ^l(x) 
described by eq.dlj) and eq.(|l5|) even for n > —1. In this 
latter case, our results are no longer exact but we can ex- 
pect that they should still provide a good approximation 
to the actual pdf. 

3.2. Generating function 

Using the spherically symmetric saddle-point obtained in 
Sect. 3.1 we can use a steepest-descent approximation to 
the path- integral (12|). That is, we approximate this path- 
integral by the Gaussian integration around the spherical 
saddle-point (5l(x). This yields: 



ni/2 



,Piy) c (DetA^i) ' (DctA/,) 



-1/2 -Sy/a\R) 



(23) 



where the matrix My is the Hessian of the exponent eval- 
uated at the saddle-point Sl obtained in Sect. 3.1: 



Mj^(xi,X2) = 



d^S/a^) 



,5(5l(xi))^(^l(x2)) 
S'{Sr.) 



a^R) (5(<5l(xi)),5((5l(x2)) 



Ai'(xi,X2) 



(24) 



In eq.(p^) we used the fact that the second derivative of 
Slj, is also the second derivative of pfj since Pr — i + 5lj,. 
The quantity 5*^, which appears in eq.(p3|) is the value of 
the action S[Sl] defined in eq^(0) at the spherical saddle- 
point 5l(x) derived in Sect. 3.1. As in paper II, after we 



substitute eq.([15|) into eq.(|13D we can write Sy as the so- 
lution of the implicit system: 



-y G'p{t) 



Sy^y Gpir) 



T 

Y 



(25) 



where we introduced the variable t and the function Gpir) 
defined by: 



T{h,RL) 



^L.Ri 



a{R) 



a[Tp{5L,R,Y/^R] 
and: 

Gp{t) =^p{^L,Rl) ^ PR. 



(26) 



(27) 



Using eq.(|26|) we see that the function Qp{t) deffired in 
eq.(|27|) also obeys the implicit equation: 



Qp{r) = ^p 



^[gp(r)^/^i^] 
a{R) 



(28) 



A s for Tp, the function Qp differs from the one introduced 
in paper II by a factor -f 1 because we study here pR rather 
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than Sff. For a power-law power-spectrum P{k) ex k''' 
have a{R) ex i?-("+3)/2 so that eq.(||) simphfies to: 



we 



Gpir)=Tp 



-r Gpir) 



-(n+3)/6 



(29) 



The previous equations allow us to evaluate the gener- 
ating function ip{y) through e q.(p3|). H owever, as for the 
quasi-linear regime studied in [paper II we must first en- 



sure that the spherical saddle-point derived in Sect. 3.1 
indeed the global minimum of the action S'[(5l]. Since we 
investigate here another regime (p^ — + and not a —f 0) 
we have to consider this point anew. Besides, we need to 
make sure that the Gaussian approximation really makes 
sense. Indeed, contrary to what occurs in the quasi-linear 
regime, the action S[Sl] is not multiplied by a large fac- 
tor in the exponent in eq.([l^) (we do not study the limit 
l/cr^ ^ oo) so that it is not obvious a priori that the path- 
integral should be dominated by a small neighbourhood 
of the saddle-point Sl in the limit p^ ^ 0. 

3.3. Validity of the steepest-descent method 

First, we note that the regime of rare voids we investigate 
here corresponds to large positive y. This is obvious from 
the definition (||) which clearly shows that for y -^ -|-oo 
we mainly probe the behaviour of the pdf V{pr) for small 
overdensities pR ^ 0, that is for very underdense regions. 
Of course, rare voids also correspond to Sl r^ ~^ — oo, 
hence to r ^ -|-oo as can be seen from eq.(26). Moreover, 
the behaviour of J[p(Sl.Rl) i^ this regime was derived in 



eq.(21). Using eq.(E9|) this yields: 



+00 : PR ^ t/p(r) 



20 

27 



-6/(l-n) 



(30) 



Then, using eq.(|25|) we obtain 



y 



20^^-3/(4-n) ^ g ^ (l-«)/(8-2n) 



27 



1-n 



(31) 



As explained above, in order to justify the steepest- 
descent method we first need to show that the spherical 
saddle-point derived in Sect. 3.1 is the global minimum of 
the action S'[(5l], that is, that there is no deeper minimum. 
For positive y the action S[Sl] written in eq.dlS) satisfies 
the lower bound: 



y > : S[Sl] > 



'-(R) 



SL.AV'.Sr > 



(32) 



since the density pR is positive. Moreover, the kernel 
A^"'^ is positive definite, as shown by eq.(p[). Hence the 
action SISl] goes to -|-cx) for large linear density fields 
|(5l| ^ -|-oo. Since we obtained only one spherical saddle- 
point in Sect.BJ (there i s only on e solution r(j/) of eq. (psf) 
for positive y, see also paper H ) we can conclude that 
this saddle-point is also the global minimum of the action 
S[Sl] restricted to spherical states Sl- Then, we need to 
show that there is no deeper minimum realized for a non- 
spherical state S^. Unfortunately, since we do not know the 



explicit form of the functional ^^[(5^] we cannot prove this 
in a rigorous fashion. Hence we shall simply assume that 
the spherical saddle-point derived in Sect.pj is indeed the 
deeper minimum of the action. Note that this assumption 
is actually quite reasonable. Indeed, for a given "magni- 
tude" of the linear density field 5l we can expect spherical 
states to be the most efficient to build very underdense 
spherical regions. 

The second point we need to check in order to val- 
idate the steepest-descent method is to show that the 
Gaussian integration we performed in eq.(|23|) is justi- 
fied. In other words, we must show that in the limit 
y —^ -|-oo (which also corresponds to pR — > 0) which we 
study here, the action S[Si] in the path-integral ( p^ ) can 
be replaced by its Taylor expansion up to second-order. 
Unfortunately, this is a rather difficult task since we do 
not know explicitly the functional pr[Sl]. Therefore, it 
is not easy to estimate the high-order derivatives of the 
action S \Sl] defined in eq.([l3[). Moreover, as noticed in 
paper H it is not obvious that one could use standard 
perturbation theory around the spherical saddle-point Sl 
to estimate these higher order derivatives. Indeed, one can 
check that using the form of the n— order term 5'^"^(x) 
obtained from the expansion of the non-linear density 
field as a power-series over the linear density field Sl 
leads to divergent quantities (see also paper V ). This is 
similar to the usual divergences one encounters in high- 
order terms ("loop corrections") obtained from such ex- 
pansions to compute the non-linear two-point correlation 
(e.g., Scoccimarro & Frieman (1996)). Therefore, we shall 
only show that the steepest-descent method is justified 
with respect to the integration over the one-dimensional 
variable t. That is, we show that it is valid for the quan- 
tity: 



V^i(2/) 



with: 



dr 



~Si{T)/a^ 



2tt(t 



Si{T)^ygpiT) + 



(33) 



(34) 



The ordinary integral ( p3D is a simplified version of 
the path-integral (|lj), where we replace the infinite- 
dimensional variable <^l(x) by the real variable r. It also 
exhibits a saddle-point Tc which is still given by eq.(P5|) 
and the value of the exponent at this point is again given 
by eq.(g5|) with S'i(tc) — Sy. For large y we can use the 
asymptotic behaviour ( pO|) and the saddle-point Tc is given 
by eq.(|3l|). Then, making the change of variable r — t^x 
we can write eq.(p3h as: 



V'i(2/) 



Ax 



-6/(1-") 1 1^2 



+ i:r=]/^= 



27rcr 



(35) 



where we restricted the integration over cc > since for 
large positive y the integral (^) is dominated by large 
positive r. Then, it is clear from eq.(^) that the steepest- 
descent method yields exact results in the limit y — > -t-oo. 



6 
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which corresponds to Tc — *■ +c». We can expect this be- 
haviour to remain vaHd (for the most part) for the path- 
integral (12). Therefore, we shall assume in the following 
that the steepest-descent method is indeed justified. 

3.4. Asymptotic form of the generating function 

We can now evaluate the generating function ^{y) ob- 
tained in eq.(23). This expression can also be written: 



V'(2/) =X'"i/2g-s„/^^(fl) with I? = Det(AL.Afj,). (36) 



Using eq.(|25|) and eq.(30) this yields the asymptotic be- 
haviour for large y: 



y 



Sy = ay 



for the minimum of the action 5*. with: 



(l-n)/(4-n) 



-6/(4-n) 



and: 



U! 



(37) 



(38) 



(39) 



Next, we need to evaluate the determinant V. However, 
since it involves the second-order derivative of the func- 
tional /Oij[(5L], see eq.(P4|), we do no t know its explicit ex- 
pression. Therefore, as in paper II (App.B) we shall use 
the analogy with the simple integral (p3|). The steepest- 
descent method yields for ipi{y) the expression: 



V'i(y) 



1 



i + yG';iT) 



O-SyjO 



(40) 



where t is the saddle-point given by eq.(Eq). Hence in 
the case of the path-integral (O) we simply have the sub- 
stitution 1 + yQ'pij) — > T). Thus, in a first step we may 
approximate the determinant T) b y the sim ple quantity 
1 + yQ"p{j)- However, as noticed in paper II this analogy 
does not take into account a physical process which is 
specific to the non-local problem of large-scale structure 
formation: the stronger expansion of underdense regions. 
Indeed, the ordinary integral (|3^) only applies to a local 
physics. On the other hand, while their density contrast 
decreases towards — 1 underdense regions also depart from 
the mean background expansion and they actually grow in 
comoving coordinates. In fact, the volume they occupy is 
larger by a factor [R/Rl)'^ = p^^ relative to the initial co- 
moving value. This increases the pdf V{pr) as well as the 
average ?/'(y) — (e~^'"*/°' ) by the same factor. Therefore, 
we approximate the generating function ^(y) by: 



y -^ +<x : %p{y) 



1 



1 



,-Sy/a\R) 



^^•(^) ^i + yg;{r) 



(41) 



Here the overdensity pR is given by eq. ( P7| ) . It is the over- 
density associated with the saddle-point r. Note that the 
factor p'^ corresponds to a non-local physics. Indeed, as 



the underdense bubble grows in comoving coordinates it 
gobbles neighbouring regions whose density becomes gov- 
erned by the properties of this initially remote underdense 
patch. 

Note that we can actually see the factor 1/ pa appear 
in the path-integral ( |l2| ) through the following discussion. 
When we translate the spherical saddle-point (5l(x) ob- 
tained in Sect. 3.1 by a vector r we do not change the term 
((5l-A^^.(5l) in the action (|3|), see eq.(||), since a transla- 
tion only yields a phase shift to the Fourier components 
5l (k) (this is related to the fact that the initial conditions 
are homogeneous). Moreover, the overdensity pr[5l\ is not 
significantly changed as long as the displaced "bubble" of 
radius R shows some overlap with the spherical cell of ra- 
dius R centered on the origin. This holds as long as r ^ _R. 
Thus, this yields a set of linear states 5[ (x) over which 
the action S[5l\ is almost degenerate. Therefore, we must 
sum up the contributions of these various states (this is 
similar to the usual case of degenerate saddle-points). If we 
discretize the comoving coordinates x by a grid of step Aa; 
we obtain a number of such states which scales as (R/Ax)^ 
(since they are translated from the origin by a length r of 
order of or smaller than R). Next, we normalize to the lo- 
cal physics where each region follows the mean expansion 
of the universe, which yields a factor (R/Rl)^ = 1/pr 
(the discrete step Ax eventually disappears as it should). 
Therefore, we naturally recover the prefactor 1/ pr. 

This mechanism also shows that we could have some 
additional deviations from the prefactor written in eq. ( [4l| ) 
if the dependence of the action on the magnitude of Sl is 
not the same along all directions (i.e. there are almost flat 
or very sharp directions). However, this appears rather 
unlikely. Indeed, as shown in eq.(^5[) the variation of the 
action around the saddle-point is governed at the same 
order (in the simple case) by the functional /o_r[(5l] and 
by the simple quadratic term (Sl-A^^.Sl). This last term 
only remains constant for translations (treated above, 
which yield the factor 1/pr) and for rotations. However, 
this last symmetry only gives a constant factor in which 
is absorbed into the normalization of the path-integral. 
Therefore, the quadratic term (Sl.AJ^^.Sl) should prevent 
any new "flat direction" . On the other hand, there might 
exist "sharp" directions if the density pR were strongly 
unstable with respect to non-spherical perturbations, for 
instance. However, this is not the case here since the dy- 
namics is actually stable around the saddle-point. Indeed, 
it is well-known that in the similar case of the expan- 
sion of low-density universes (i.e. with fim < 1) density 
perturbations do not grow (i.e. the linear growing mode 



-D+(t) saturates as soon as ^m ^ 0.2), see Peebles (1980) 



However, we can clearly expect that an exact calculation 
would give a multiplicative numerical factor of order unity 
with respect to eq.(^l|). In any case, note that the expo- 
nential term is exact (for —3 < n < —1) since it only 
depends on the value of the action at the saddle-point. 
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3.5. Low-density tail of the pdfV{pR) 

From the generating function "ipiy) obtained in eq.([4l|) 
we can now derive the low-density tail of the pdf V{pr). 
Indeed, from the inverse Laplace transform (W and the 
definition dll]) we can write: 



-PipR) 



-\-ioo 



dy 



^ypnl'y'^{R) 



Using eq.(|4l|) this yields: 



V'(2/)- 



(42) 



-PipR) 



AyPa~Sy]/a 



-- 2^ja2 ^^(^) 



1 



yG';{r) 



(43) 



Then, in the limit p^ — > we can evaluate this integral 
by an ordinary steepest-descent method which gives: 



Vipi- 



1 



1 



.-^V(2t^) 



27rcr 



PB. 



^/i + y5';iT) \F^ 



(44) 



Here the variable t is again given by Gp{t) = pR while Sy 
is the second-derivative with respect to y of the value of 
the action 5*1, at the saddle-point associated with r. From 



eq.(25) we get: 



dr 



Q»^y^ l + vQ;{r) = -Q'p{r\ 



Sy 



SO that eg. (|44|) also writes: 
1 1 1 



dy 



V{pr) 



-V(2-=) 



^(T PR \Q'p{r)\ 
Then, using eq.(p3C|) we obtain 
1 1 - n 27 



(45) 



(46) 



V{pi- 



V2^cr 6 20 



-13 

Pr" e 



'^ p-(i)^P«'^-"'^V(2^=) 



(47) 



The asymptotic expressions (^^ and (^) hold for very 
rare underdensities, beyond the cutoff r^ or p^ of the pdf 
T^iPR) ("v" for voids). This characteristic underdensity 
is set by the condition r ^ o-{R), see eq.(46). In the 
quasi-lin ear regime we have r ~ ~Sr, see eq.(P6|) and 
paper II, which gives 6y ^ —a{R) as expected. In the non- 
linear regime where u{R) 3> 1 we have r„ ^ 1 so that 
we can use the asymptotic behaviour (pQ) which yields 
Pt, ~ cr(i?)~^/(^""'. Therefore, the typical overdensity p^ 
of voids is given by: 



a{R) 



-6/(l-«) 



(48) 



o- <C 1 : Pu = 1 — cr(i?), cr > 1 : /9„ 
In the non- linear regime this yields: 
a{R) > 1 : p„ == (7(i?)-6/(i-") ex i?3(«+3)/(i-«)_ (49) 

Thus, the expressions (E6|) and ( ^7| ) provide the asymp- 
totic behaviour of the low-density tail of the pdf VJpr). 
They apply to pR ^ py. From eq.(pO|) and eq.(31) we 
see that the density py corresponds to the value y^ of the 
Laplace variable y, with: 



a :S> I : yv ^ Pv^"^' 



■0/3 ^ ^(8-2«)/(l-n)^ 



(50) 



In particular, the asymptotic form ( p4| ) of the generating 
function ijj{y) actually applies to y 3> 2/t,. 

Here, we must stress that the results ([4q ) and (47) 
have been directly derived from the equations of motion. 
For n < — 1 the exponential in eq.(p7|) is exact but the 
prefactor is only approximate because of the approxima- 
tion (^) we used for the determinant V. For n > — 1 the 
normalization factor in the exponential is also approxi- 
mate because we did not use the exact saddle-point, as dis- 
However, the characteristic exponent to 



cussed in Sect. 3.1 



-(l-«)/3 



in eq.(|37|) and the power pp, within the exponential 

cutoff in cq.(|47|) should still be exact. Besides, we recall 
that cq.(|47|) holds for arbitrary values of the variance a, 
provided pr-^ Pv Indeed, we only used the limit of very 
underdense regions (rare voids) and the actual value of a 
was irrelevant in the calculation. In fact, the normalization 
of the power-spectrum does not even appear in the char- 



acteristic action S[5l\ defined in eq.(13) which describes 
the physics of our system ! Thus, our result (E^ is valid 
from the quasi-linear regime up to the highly non-linear 
regime. The influence of a only shows up in the constraint 
PR ^ Pvi as shown by eq.(p8|). Indeed, it is clear that 
as time goes on and gravitational clustering builds up the 
characteristic overdensity pR of typical voids evolves. That 
is, the cutoff at low densities of the pdf 'P{pr) is steadily 
pushed towards lower densities as gravitational clustering 
proceeds. 

The quasi-linear regime was already studied on its own 
in paper 11. In that previous work, we investigated the 
limit (7^0 which gave the pdf V{pr) (or V{5r)) for all 
density contrasts 5r, in the quasi-linear regime. Indeed, in 
this limit any finite density contrast becomes a rare event 
as soon as a <C \5r\. Since the calculation involved the 
same spherical saddle-point ( p5| ) described in Sect. 3.1 the 
expression (E^) is consistent with the results obtained in 
that previous paper. Of course, the main interest of eq.(E^ 
is that it also applies to the non- linear regime a > 1. 
Indeed, this regime is more difficult to handle and very 
few rigorous results were known so far. Moreover, in order 
to build significant underdensities one must be at least in 
the mildly non- linear regime cr ^ 1. 

3.6. Comparison with previous works 
3.6.1. Perturbative methods 

We can note that the approach described in the previ- 
ous sections bears some similarity with a perturbative 



study developed in Bernardeau (1994). This author inves- 
tigated the mean non-linear dynamics in the limit of rare 
events (i.e. large density fluctuations) through a pertur- 
bative method which assumes that the non-linear density 
field can be written as a perturbative expansion over the 
linear density field. Note that this problem is not obvi- 
ous a priori. Indeed, although the mean density profile, 
with the constraint that the average linear density con- 
trast within the radius Rl is some fixed value Sl.r^ : obeys 
eq.dlq) (i.e. exactly the profile of our saddle-point), it 
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does not follow that the mean non-linear profile should be 
given by the non-linear evolution of the mean linear pro- 
file (i.e. these two operations may not commute). However, 
Bernardeau (1994) noticed from a perturbative treatment 
that this is actually the case as one eventually recovers 
the usual spherical dynamics. This agrees with our results 
(Eq) and (p7|). In fact, our approach provides a simple ex- 
planation for this feature. This peculiar spherical solution 
of the dynamics is actually a saddle-point of the action 
S[6l] so that it governs the tails of the pdf V{pr). Note 
also that our method is much more intuitive and simpler 
as it clearly reveals the underlying physics. 

Next, we stress that our method is non-perturbative 
and it yields exact results (as long as one can identify 
the exact minimum of the action). By contrast, perturba- 
tive methods are based on the expansion of the non-linear 
density field over the growing linear density field, which 
is then plugged into the hydrodynamical approximation 
of the equations of motion. Therefore, such approaches do 
not provide complete proofs since the perturbative expan- 
sions should diverge. Moreover, they ca nnot go beyond 
shell-crossing. Thus, Bernardeau (1994) noticed that for 
power-spectra with n > — 1 the first correction (i.e. sub- 
leading term) obtained by the perturbative approach di- 
verges. This implies that the perturbative method fails 
for n > — 1. This feature can easily be understood from 



the discussion given in Sect. 3.1. Indeed, we noticed that 
for n > — 1 the spherical saddle-point of the action ex- 
periences some shell-crossing for i?^ ^ R^. Therefore, it 
cannot be obtained by perturbative means and all pertur- 
bative methods must diverge. Nevertheless, this does not 
invalidate our steepest-descent method which is essentially 
non-perturbative. It merely means that it is more difficult 
to obtain an analytical expression for the exact m inim um 
of the action S[5l\- Besides, as described in Sect.^ our 
approach provides a very convenient way to obtain ap- 
proximate results in such cases. We simply need to use an 
approximation for the spherical saddle-point. As argued 
in Sect.SJ. we can expect this procedure to provide very 
good results for the case of rare underdensities since even 
for n > — 1 shell-crossing only involves a limited range of 
radii along the density profile of the saddle-point. In par- 
ticular, we expect that all characteristic exponents (e.g. 
the power pj^ in the exponential term in eg. (p7|) ) 

should remain exact. 

Note that this feature for —1 < n < 1 clearly shows 
once more that perturbative results should be viewed 
with caution. Indeed, at leading-order the perturbative ap- 
proach yields the usual spherical dynamics, as in eg. ([l7|) . 
However, as we explained in Sect.^AJ this does not cor- 
respond to the exact spherical saddle-point which means 
that the leading-order behaviour of the pdf Vipn), or of 
the generating function '0(2/) j is not given by this sim- 
ple expression. Indeed, we can expect the actual min- 
imum of the action S[Sl] to be slightly different from 
the value computed from eg.(|l4[), which translates into 
a slightly different value for the numerical factor a in 



eg.(|37|). Therefore, the divergence of the subleading terms 
actually leads to a correction to the finite leading term 
derived from the perturbative method. Then, we note 
that for hierarchical scenarios all perturbative series ac- 
tually diverge because on small scales the density field 
is always a non-perturbative guantity (see paper 1 and 
paper V). Moreover, it is well-known that in standard 
perturbative expansions one actually encounters diver- 
gent guantities beyond some finite order ("loop correc- 
tions", e.g. Scoccimarro & Frieman (1996)). As a conse- 
guence, there can be no guarantee that perturbative re- 
sults (even though finite and restricted to leading-order 
terms) make sense. The only way to obtain firm results is 
to use non-perturbative methods which can overcome (at 
least in pr inciple) t hese problems. The goal of our previ- 
ous work ( paper H ) and of the present article is precisely 
to develop such tools. 

Finally, since we shall investigate the high-density tail 
of the pdf in Sectjijwe note h ere that perturbative meth- 
ods as in Bernardeau (1994) are restricted to the early 
non-linear stages of the dynamics before shell-crossing oc- 
curs. By contrast, the high overdensities we shall study 
correspond to non-linear objects which have already viri- 
alized and where shell-crossing plays a key role. 

3.6.2. Spherical model 



As in [paper II| , we note that the expression ( [IGj) can ac- 
tually be recovered from a very simple spherical model, 
detailed for instance in Valageas (1998). This phenomeno- 



logical model rests on the approximation: 



dS{l + 6)ViS) 



dSLVLiSL). 



(51) 



Here Vl{5l) is the linear pdf of the linear density contrast 
5l within a spherical cell. This relation merely states that 
the fraction of matter contained within spherical cells of 
radius R and non-linear density contrast larger than 5ji 
is approximately egual to the fraction of matter which is 
enclosed within spherical cells of Lagrangian radius Rl 
and linear density contrast greater than Sl,Rl- Here Rl 
and Sl,Rl ^re related to R and Sr by the usual spherical 
dynamics, as in eg.([17|). Note that this is very close to 
the Press-Scheehter prescr iption without the factor 2, see 
Press fc Schechter (1974) , Next, we note that the linear 



pdf VLi^L) at scale R exhibits a simple scaling over the 
variable v through: 

5l 



Vl{5l) d5L = V^^\^) d;^ with 



a{R) 



and: 



^-u-/2 



'271 



(52) 



(53) 



Substituting eg.(M) into eg.(| 
respect to 5r we obtain: 

^ ,r N 1 I dv 

Vs{Sr) 



and differentiating with 



1.^2 



1 + Sr 727? d6R 



(54) 
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with: 



OL.Rl 



a{Ry 



(55) 



Here the subscript "s" refers to the "spherical" model. 



eq.(55) we introduced the variable t defined in cq.(& 
Moreover, using the function Qp{t) introduced in eq.(t 
we have: 






1 



dr 



1 



1 



<j{R) A5r a{R) mT)\ 



(56) 



Then, substituting cq.(M) into eq.(p4) we recover the ex- 
pression (|4^). Thus, for low densities pR <C Pv the spher- 
ical dynamics correctly describes the leading order be- 
haviour of the pdf V{pr). This is not surprising: it is sim- 
ply due to the fact that the saddle-point of the action 
S[5l\ obtained in Sect.3T is spherically symmetric. This 
holds because the initial conditions are homogeneous and 
isotropic and we study the statistics of the mean overden- 
sity PR over a spherical cell of radius i?, centered on the 
origin (for instance), which preserves the spherical sym- 
metry of the problem. Note however that the determinant 
T> which appears in the prefactor of the generating func- 
tion V'(2/) in eq.(p6|) takes into account the deviations of 
the initial conditions from spherical symmetry (at leading 
order in the limit p^ -^ 0). As discussed in Sect.^.4[, we 
can expect the exact Gaussian integration over the non- 
spherical density fluctuations 5l around the saddle-point 
to give a multiplicative numerical factor of order unity 
with respect to eq.(p^). Moreover, as shown in Sect. ^.l| 
for — 1 < n < 1 the spherical model only yields an ap- 
proximation for the exponential cutoff of the low-density 
tail since the saddle-point given by eq.(^ is no longer 
exact. 

As already noticed in Valageas (1998)| , the overden- 
sity py obtained in eq. (p^) corresponds to density fluctu- 
ations which would occupy all the volume of the universe 
at the time of interest. This is the typical density of voids 
which fill almost all the volume of the universe in the 
non-linear regime. It is clear that for higher densities the 
deviations from the spherical dynamics and the effects of 
shell-crossing play a key-role and they must be taken into 
account. 

3.6.3. Non-linear scaling model 

Here, we point out that our results ( |37| ) and (|4^) are remi- 
niscent of the prediction of a non- linear hierarchical ansatz 
investigated in Balian fc Schaeffcr (1989) . This model is 
based on the assumption that in the highly non-linear 
regime (cr 3> 1) the non-linear many-body correlation 
functions Cg(xi, ..,Xg;a) obey the scaling law: 

C,(Axi,..,Ax,;a) =a3(9-i) A^^^'^-D 4(xi,..,x,) (57) 

for arbitrary A > and any time. Here a{t) is the scale- 
factor while 7 is the slope of the non-linear two-point cor- 
relation function ^. This scaling law can be derived from 



the stable-clustering assumption (Peebles (1980)). In this 
case, for a power-law power-spectrum P{k) ex fc" we have: 



3(3 -hn) 



7 = 



Then, it is convenient to introduce the quantities: 



(58) 



I, 



C 



- with ^R) 



dxi..dXg 



Cg(xi,..,Xq) (59) 



where we note ^ = ^2- Thus, the parameters Sq yield the 
cumulants (/9^)c = S,„ (with ^^ = 1). Besides, the scaling 
laws (^^ imply that the coefficients Sq do not depend on 
scale nor on time. Next, the Laplace transform ip{y) of 
the pdf defined in eq. (||) is related to these cumulants by 
the standard property (see any textbook on probability 
theory): 



^V;(y)]=^^^(p|,)e2/'^. 
Using eq.(|59|) this yields: 



V'(y) = '0(yO with tpiy) = e 
and: 

q=l ^■ 



^ p-viv)/i 



(60) 



(61) 



(62) 



Note that within this model the generating function (p{y) 
is scale and time independent. Besides, it provides the pdf 
VipR) through eq.(||) which reads: 



+ioc 



r{pR) = I -^ elyPH-v{y)]/i 

27ri^ 



(63) 



where the tilde "~" refers to the hierarchical ansatz. As 



argued in Balian fc Schaeffer (1989) , it is natural to ex- 
pect a power-law asymptotic behaviour at large y: 

y -^ +00 : (p{y) ~ a y^""^ with a > 0, < lu < 1. (64) 



Note that this is similar to our results ( |36| ) and (pTh for 
the generating function V'(2/)- Apart for the factor V in 
eq.(pq), the only difference is that within this hierarchi- 
cal ansatz the linear variance a^ must be replaced by its 
non- linear counterpart ^. From eq.(|63|) and eq. (|64|) one 
obtains for small overdens ities pR <^ ^ the behaviour 
( iBahan fc Schaeffer (1989)|) : 



Pi?<?: V{pR.)=a '^- > £, gc{z) 

where we introduced the variable z defined by: 

z = a '"^ > E, PR 

and the function gcj{z) given by: 



gc^iz) 



-\-ioc 



dt 

27ri 



(65) 



(66) 



(67) 
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The function gcj{z) exhibits a sharp cutoff for 
small z which can be computed by an ordinary 
steepest-descent method. This eventually yields 



(Balian fc Schaeffer (1989)) 



z « 1 : V{pr) 



a-i/(i- 



T 



cO/(1-lD) /(l-(D)i/'^ 



2ttuj 



X 2-(l+i')/(2c:') g-i;[z/(l-tD)]-'i- 



.(68) 



At first sight, this cutoff looks similar to eg. ([47|) . First, 
we note that we can make the typical void density p^ 
implied by eg. (pSJ) to coincide with our result (^. Indeed, 
the value of the parameters a and uj is not predicted by 
the hierarchical ansatz ( p7| ) . The low-density cutoff of the 
pdf V{pr) is given by z ~ 1, which yields for the typical 
overdensity p^ of voids: 



a(i?)»l: p„=r^^^'"^^oci?^^/(i^ 



^) 



Using eg.(pq), the comparison with our result 
gives: 



OJ 



6 



(69) 
for py 

(70) 



Here, it is interesting to note that the value ([7C|) for ui 



was already obtained in Valageas & Schaeffer (1997) from 



the stable-clustering ansatz coupled with the spherical col- 
lapse model. Let us recall brieffy the main properties of 
this model. As in Sect. 3.6.2, it is based on the approx- 



imation (|5l| ) which is used to relate the non-linear den- 
sity field (5(x) to its linear counterpart Sl (x) . Then, in 



Valageas fc Schaeffer (1997) we used eg.(|5l|) to "derive" 
the pdf V{pr) in the non- linear regime. In particular, we 
considered collapsed objects which have already virialized. 
Then, within the framework of the stable-clustering ansatz 
we have pr cx p~^ (x a'^ while the linear density contrast 
obeys (Sl(M) cx a (in a critical density universe) which 
yields pR cx S^ and ^{Sl) oc 5£. Next, substituting this 
behaviour into eg.(Bl|) we obtain: 



Pv < p/? < C 



ViPR) (X 32 ( ^ 



(71) 



where the exponent w is given by eg. (|70|) and the low- 
density cutoff is given by eg. (^9|) . Note that the power-law 
behaviour (jnh agrees with the scaling (^) over the range 
pv ^ PR ^ S, where both formulae overlap. Indeed, for 
large values of z the functio n gcbjz) defined in eg.(|67|) ob eys 
the asymptotic behaviour ( Balian fc Schaeffer (1989) ): 



z'>l ■■ gcj{z) 



1 



nQ) 



(72) 



This is obtained by expanding the term —t^^'^ in the ex- 
ponent in eg.(|67|) since for z :^ 1 only small values of t 
contribute to the integral. Here, we must point out that 
this line of reasoning involves virialized objects which sat- 
isfy pR ^ yO«, whence the lower bound for the range 
of validity of eg. ([Tl]) . In particular, at this stage stable- 
clustering does not imply the scaling (|65|) for pR ^ p^. 



In order to obtain eg.(B5[) down to pR -^ one needs to 
assume that the power-law behaviour (B3) extends up to 
y —t +00 at all times and scales (the function (p{y) is scale 
and time- independent in this framework), as assumed in 



Balian & Schaeffer (1989). However, this is clearly incon- 
sistent with our rigorous results (Eq) and (|47|). Indeed, 
using eg.(^fl) we obtain (1 — a))/J; = (1 — n)/(5 -I- n). 
This implies from eg.(p8|) that the pdf V{pr) exhibits a 
low-density tail of the form V{pr) ^ e"'^''^/^"-* 
which disagrees with eg.([47|). Thus, our study explicitly 
shows that the non-linear hierarchical ansatz (m) does not 



describe the pdf P{pr) for rare underdensities pR <C Pv 
In other words, the very low density tail of the pdf cannot 
be derived from the stable-clustering ansatz. 

Actually, this is not surprising in view of the physics 
which lies behind the derivation performed in Sect. 3.1 and 



Sect. 3.2. Indeed, we have shown that the low density tail 
of the pdf (i.e. pR <^ py) is governed by the dynamics 
of spherical very rare "low-density bubbles" which are 
still expanding. In particular, in this asymptotic regime 
the expansion of the outer shells is almost "free" , that 
is their physical radius grows as i? cx i which implies 
PR (X t^^ (X 6^ ' (if rim = 1) in agreement with eg.(pi|). 
Indeed, the gravitational pull from the inner regions be- 
comes negligible. Note also that virialization processes do 
not show up at all. Therefore, we could have expected the 
stable-clustering ansatz (^7|) to be irrelevant to the be- 
haviour of the pdf V{pr) in this very low density regime. 
Nevertheless, the stable-clustering ansatz can be made 
to recover the correct void density py in eg.(^) because 
within this framework these underdense regions have just 
stopped their expansion and they "virialize" at the time 
of interest. Hence the stable-clustering assumption does 
not have any influence on their properties yet. 

Note however that this result does not imply that the 
scaling laws (|5^) are wrong. Indeed, it is clear that the 
many body correlation functions ^g involved in eg.(p7[) 
are dominated by the high-density regions pR <; ^, which 
may be governed by viriahzation processes. Our result (47) 
merely means that the low-density tail of the pdf (i.e. 
PR '^ Pv) depends on the detailed behaviour of the many 
body correlation functions ^q which is not fully captured 
by the lowest-order asymptotic behaviour (|57|). Indeed, 
the discrepancy between eg.(p^ and eg.(|65|) means that 
the power-law behaviour ( |64|) only applies up to y ^ y„ at 
most, with: 



Vv = Pv^^^ 



(73) 



as can be obtained from eg.(p3). In the highly non-linear 
regime we have pv ^ hence we get y^ ^ oo. Therefore, 
it is clear that the behaviour of <f{y) for y > jjv cannot 
be described by the asymptotic scaling laws (|5^), even if 
the latter are valid, since in the highly non-linear limit 
(7 — > -t-cx), which defines the limiting generating function 
(p{y) written in eg.(|62|), this regime disappears as j/i, is 
repelled to +00. 
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3.6.4. Numerical simulations 



Finally, our result (|£^) should be compared with numer- 
ical simulations. However, this is rather difficult since 
the steepest-descent method only applies to the far low- 
density tail of the pdf: pji <C Pv In practice, in numeri- 
cal simulations one mainly observes a sharp cutoff below 
some characteristic overdensity pn and it is not easy to 
measure with a good accuracy the shape of the pdf be- 
yond this cutoff. In fact, because of discrete effects (i.e. 
the limited number of particles) one does not really probe 
the low-density tail described by eq. ( [47| ) in current numer- 
ical simulations. Indeed, let us note V{N) the probability 
to have N particles within a spherical cell of radius R in 
a given simulation. If we assume (as is usually done) that 
discretization only adds a Poisson noise the pdf V{N) is 
obtained from its continuous counterpart V{pr) by: 



V{N) 



<ipR. VipR, 



{prN) 



N 



-PrN 



The kernel in the integrand in eq. ( |74| ) is simply a Poisson 
law and we note N the mean number of points within a cell 
of radius R. Next, we note Vo = V{N ~ 0) the probability 
to find an empty cell in the simulation. Then, from eq.([74[) 
and the inverse Laplace transform (Ba) a straightforward 
integration over pn and next over y yields: 



Vo = ^j{Na^). 



(75) 



Therefore, we see that Vq probes the generating function 
ip{y) at yo — Na'^. Moreover, it is clear that the numerical 
simulation cannot probe the continuous generating func- 
tion ip{y) to larger y, that is the pdf V{pr) to smaller 
densities. Thus, the simulation only tests the low-density 
tail ( [47| ) if this value j/o is much larger than the value yv 
obtained in eq.(pQ) which marks the onset of the regime 
associated with these very underdense regions, in the non- 
linear regime cr ^ 1. Using eq.(pfl) this constraint reads: 



a>l: yo>2/. if iV>(7'^/(i-"). 



To check whether this condition is realized in practice we 
can look at the typical numbers reached in current sim- 
ulations. For instance, the numerical simulations used in 
Valageas et al. (2000)| contain 128^ ~ 2 x 10^ particles 
within a box of size 256 Mpc (the length scale is actually 
arbitrary for a scale- free power-spectrum). At the end of 
the simulation (when the largest scales approach the non- 
linear regime) the value cr = 10 (not to take a too large 
number for the r.h.s. in eq.([76[)) has only been probed in 
a reliable way by the scales i? ^ 4 Mpc (the scale 8 Mpc 
may have reached ct = 10 but it is not accurate because of 
finite size effects). The mean number of particles within 
a cell of radius 4 Mpc is A^4 ~ 8.4 while for n — —2 (the 
most favourable case) we have 10^/*^^"") — 100. Therefore, 
the constraint (|7q ) is very far from being satisfied. Thus, 
we can conclude that the low-density cutoff seen in cur- 
rent numerical simulations is governed by discrete effects 
and it does not probe the actual low-density tail ([47|) of 



the continuous pdf V{pr). This requires a much larger 
number of particles. 

Finally, we note that numerical simulations have been 
used to check the scaling model (p^, see for instance 
polombi et al. (1996) . To do so, one uses the fact that 
within this framework the probability Vq to find an empty 



cell can be written (e.g., Balian & Schaeffer (1989)): 



Vo 



-V(WC)/C 



(77) 



This relation can be obtained in a manner similar to the 
derivation of eq. ( [75| ) . This allows one to derive the expo- 
nent a) defined in eq.(p3). As explained above, our results 
do not confirm nor invalidate these works since the regime 
probed in these simulations is not covered by our present 
study: it corresponds to "high densities" pn ^ py above 
the range of validity of eq. (E^ . 



(74) 4. Rare overdensities 



Finally, we investigate in this section the high-density 
tail of the pdf V{pr). This problem is more difficult 
than the study of rare voids since shell-crossing now 
plays a key role. In particular, we do not know the ex- 
plicit analytic expression of the functional pr[Sl], even 
when we restrict ourselves to spherically symmetric states. 
Therefore, we did not manage to derive the asymptotic 
behaviour of 'P{pr) for large pR in a fully rigorous man- 
ner. Nevertheless, we shall discuss the expected proper- 
ties of this high-density tail in the spirit of the steepest- 
descent method used in Sect.pl Here, by rare overdensities 
we mean massive halos which have already collapsed, that 
is we consider highly non-linear objects where the local dy- 
namical time (over the radius R) is smaller than the age 
of the universe. Thus, for approximately spherical halos 
the particles have already undergone several oscillations 
through the cluster. 



(76) 4.1. Spherical collapse 



Since we look for a spherical saddle-point of the action 
S[5l\ written in eq.([l3|) we first recall in this section the 
non- linear dynamics of spherical states. More precisely, we 
consider linear density profiles of the form: 



Sl.r oc Af ~' (X R 



-3e 



(78) 



This problem was studied in Fillmore fc Goldreich (1984)| 
and we briefly recall below their main results. The advan- 
tage of such power-law linear density profiles is that the 
evolution is self-similar. Indeed, there is no characteristic 
length scale or time in this problem so that the dynamics is 
self-similar. Thus, the system at a later time is equivalent 
to the same system seen at a smaller scale: a rescaling in 
time can be absorbed by rescaling the length scales. This 
allows one to explicitly solve the dynamics. For in stance, 
the case e — 1 was studied in Bertschinger (1985) from a 
Lagrangian point of view, where one follows the trajectory 
of individual particles (or spherical shells). This trajectory 
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r(ri,t) depends on time t and on the initial radius of the 
particle at some fixed initial time ti (or equivalently on 
the mass located within this spherical shell in the linear 
regime). Then, the self-similarity of the dynamics implies 
that all particles follow the same trajectory rescaled in 
proper units (e.g., the radius and the time of the first turn- 
around). Thus, the system is fully described by a function 
of only one variable, e.g. the time-dependence of the ra- 
dius of a particle with an arbitrary initial radius, since 
the trajectories of all other particles can be obtained from 
this one by a simple time and length rescaling. Then, this 
function is seen to obey an ordinary integro-differential 
equation which yields all properties of the system. 

The behaviour of the mass shells can be described as 
follows. After an initial stage of expansion (when 5l,Rl ^ 
1) the particle turns around at a radius rta at time ita- 
Next, the particle oscillates through the center of sym- 
metry of the system. As time goes on the mass which 
has already turned-around increases so that the parti- 
cle is buried ever more deeply within the collapsed halo. 
Besides, the particles enclosed within the central regions 
arise from a greater number of mass shells. Therefore, the 
density within these inner regions grows while the am- 
plitude of the radial oscillations of a given particle de- 
clines. However, two behaviours can occur, depending on 
the initial slope e. For sharp linear density profiles with 
e > 2/3 the amplitude of the particle oscillations asymp- 
totically stabilizes to a finite radius (which scales as its 
first turn-around radius) so that the actual density pro- 
file /9(r, t) becomes time-independent (in physical coordi- 
nates r). This is due to the fact that the contribution of 
outer mass shells to the mass enclosed within a fixed phys- 
ical radius R is negligible. Indeed, although an increasing 
number of shells "visit" this central region as they pass 
through the origin the time they spend within r < R be- 
comes steadily smaller (they spend most of their time at 
radii r of or der of their large turn-around radius). Thus, 
as shown in Fillmore fc Goldreich (1984) the asymptotic 
non- linear density profile is: 



(79) 



e>-: p(r, i) ex r-9^/(i+3^) 
o 

which does not depend on time. On the other hand, for 

shallower linear density profiles with e < 2/3 the influence 

of the outer mass shells is no longer negligible. Then, the 

amplitude of the particle oscillations no longer stabilizes: 

it slowly decreases as the mass which has already collapsed 

grows. This gives rise to a density profile which undergoes 

an adiabatic increase: 



e < o ■ Pi^^ ^) °^ ^ 



(4-6e)/(9e) ^-2^ 



(80) 



Note that the radial slope no longer varies with e but the 
time-dependence explicitly depends on e. 

4.2. Candidate for a spherical saddle-point 

Now, we can apply the results recalled in the previous 
section to the formation of rare massive halos. Since our 



system is spherically symmetric we can look for a spher- 
ical saddle-point of the action S'[(5l]. Indeed, note that a 
spherical state 5l,c(x) which is a minimmn of the action 
restricted to spherically symmetric linear density fields is 
automatically a saddle-point with respect to transverse 
directions. Indeed, assunnng the functional pr[Sl] can be 
expanded as a Taylor series around this spherical point 
Sl.c the variation of pji at linear order over a small per- 
turbation Sl = Slc + ^Sl is given by: 



ApR = / d: 



S{Sr) 



6i5Li^)) 



A5l(x). 



(81) 



Because of spherical symmetry, the first-order derivative 
6{Sr)/5{6l{^)) at the point Sl.c only depends on |x|. 
Therefore, the integration over angles in eq.(|8l|) vanishes 
for any deviation of the form x{x)Yi"^ {6 , cj)) with / 7^ 0. 
In a similar fashion, the linear deviation (J^ ^.A^^-A^^) 
which arises from the second term in the action dl3) is also 
zero. Therefore, the action S[Sl] only shows a quadratic 
variation over ASl for non-spherical perturbations, which 
implies that Sl.c is also a saddle-point with respect to 
these non-spherical directions. 

Thus, we only need consider spherical linear den- 
sity fields (assuming there are no deeper non-spherical 
minima). However, even the restriction of the functional 
Pr[Sl] to such spherical states Sl is unknown since eq.([l7|) 
breaks down because of shell-crossing. Nevertheless, in 
some cases we may still approximate the spherical func- 
tional pr[Sl.r"] (here R" is a dummy variable) by eq.(p7|). 
More precisely, this should provide a good approximation 
as long as the slope of the density profile at large radii 
i?^ > Rl is sufficiently large, that is Sl^r' oc i?^" with 
a > 2, see eq.(|7^) and eq. (|79|) . Here, the scale Rl is the 
Lagrangian scale (i.e. mass scale) of the particles with a 
turn-around radius equal to R. Indeed, for such a steep 
profile we know that the mass within the radius R is (up 
to a factor of order unity) the mass which was enclosed 
with in t his shell in the linear regime, as we recalled in 
Sect. 4.1. This justifies the use of eq.([l7|). As we recalled 
in Sect.O this yields a spherical saddle-point which is flat 
within Rl and which exhibi ts a pow er-law decline at large 
scales of the form (^9|), see [paper II for a detailed deriva- 
tion. Then, the comparison with the constraint in eq.([79|) 
implies n > — 1. 

Thus, we see that for n < —1 the linear density profile 
is too shallow to stabilize the turn-around radius of the in- 
ner mass shells. Then, the mass enclosed within any radius 
in the halo is governed by the outer shells. The non-linear 
density profile adjusts to p{r) ex r~^ so that the mean 
density within the radius R only depends on the scale Rl 
which has just turned around through pn = Pc{R/Rl)~'^, 
where pc is a simple normalization factor of order unity. 
In order to derive the spherical saddle-point in this case 
we can proceed as follows. Since in this regime the func- 
tional pli [Sl] only depends on the scale Rl where the mean 
density contrast reaches the value Sc of order unity (i.e. 
the largest non-linear scale) we first minimize the action 
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S[5l] at fixed Rl and finally wc minimize over R^. Thus, 
we must minimize S[6l] with the constraint: 



Vl 



dx 



Sl (x) = Sc 



(82) 



As usual, this is done through the introduction of a 
Lagrange multiplier A. Therefore, we need to minimize 
the action S\ given by: 

2 

Sx[6l,Rl] = ypB.iRL) + ^SL.Al\5L 



+A(^/"dx<5i(x)- 



{x < Rl) 
Vl 



- Sc (83) 



where 9(x < Rl) is a top- hat with obvious notations. 
Minimizing S\ with respect to i5l(x) yields: 



^l(x) = - 



A 



Vl 



dx^ 



Al(x,x'). 



(84) 



Note that the linear density profile is exactly of the form 
(ttq). This is not surprising since in the regime relevant 
to Sect. 3.1 the density pR does not depend either on the 
details of the density profile but only on the variables Rl 
and 5l,Rl- Hence the spherical saddle-point we obtain for 
n < — 1 is flat within the scale Rl and it decreases at 
larger scales where it is in the linear regime. 

Thus, we see that these considerations suggest two very 
different behaviours. For steep power-spectra with n > —1 
the high-density tail of the pdf V{Sji) at scale R would be 
related to large density fluctuations over the Lagrangian 
scale Rl associated with the Eulerian scale R, embed- 
ded within larger halos. By contrast, for shallow power- 
spectra with n < —1 the density fluctuations at scale R 
would be governed by the collapse of much larger scales 
Rl which are just turning non-linear. In the non-linear 
regime cr(i?) ^ 1 this scale Rl would be much larger 
than R. However, we note that numerical simulations do 
not show such a transition at n = — 1. In particular, they 
are roughly consistent with the stable-clustering ansatz 
(e.g., [Valageas et al. (2000)| , [Colombi et al. (1996)[ ) which 
would be strongly violated in case the pdf would be gov- 
erned by the saddle-point (p4|). The reason behind this 
apparent discrepancy is that the pdf V{pr) is not domi- 
nated by this spherically symmetric saddle-point. Indeed, 
this spherical state Sl only governs the pdf if the action 
remains close to this minimum over a sufficiently large 
region of phase-space, i.e. for linear density fields which 
show some slight deviations from this spherical state. For 
instance, in order to apply the steepest-descent method to 



rare underdensities in Sect. 3.2 we had to check that the 
path-integral 1l3) is really dominated by the Gaussian in- 
tegration around the saddle-point, see the discussion in 
Sect. 3^. We still need to investigate this point for high 
overdensities. As we explain below, it happens that such 
a study reveals that the path-integral is not governed by 
this spherical saddle-point. 



4.3. Virialization processes. Radial-orbit instability 

In order to check whether the steepest-descent method 
around the spherical state obtained in the previous sec- 
tion is valid, wc must study the behaviour of the action, 
hence of the functional pr[Sl], around this point. We may 
first start by investigating linear perturbation theory. As 
described below, this will prove sufficient to invalidate the 
steepest-descent method. This will also provide some use- 
ful information about the structure of collapsed halos and 
virialization processes. 

The spherically symmetric saddle-points obtained 
in Sect, 4.2 exhibit purely radial motions. Therefore, 
they give rise to non-linear spherical halos with ex- 
actly radial orbits. Such orbits are known to be un- 



stable (see for instance Palmer & Papaloizou (1987) and 



Polyachcnko (1992) for a presentation of the "radial- 



orbit" instability in some specific cases) hence we can sus- 
pect these halos to be unstable. This implies that the ac- 
tion S[5l] would increase very fast for small non-spherical 
deviations A6l around the saddle-point. Then, the den- 
sity fluctuations at scale R may not be governed by these 
exactly radial solutions of the dynamics because suffi- 
ciently spherical states are too rare. As discussed below 
this is indeed what occurs in our case. Thus, we study 
in AppJAJ the linear perturbation theory around spher- 
ical halos with nearly radial orbits. Since the growth 
rates uj we shall obtain are much larger than the typ- 
ical frequency flo ^ I /to where to is the dynamical 
time (which is of the order of or smaller than the Hubble 
time tn) the growth of the halo over the time scale ijy 
is irrelevant. Therefore, we consider static spherical ha- 
los with radial orbits. Note that this is rather different 
from the systems investigated in previous works (e.g.. 
Palmer & Papaloizou (1987), Polyachcnko (1992)| ) where 
radial orbits only involved a small fraction of the matter 
content of the halo. In particular, as explained in App.H 
while for such systems the authors found slow growth rates 
Lj ^ f2o as the perturbations develop through a resonance 
2 : 1, in our case the instability is much more violent be- 
cause it involves the whole halo and it leads to very high 
growth rates to ^ flo. 

Indeed, considering a halo with nearly radial orbits, 
or more precisely where the typical angular momentum 
p of the particles is very small, we show in AppjAJ that 
non-spherical perturbations are strongly unstable with a 
growth rate of order: 



flo \ — if M< -^0, 
p 



(85) 



where fto is the typical orbital frequency of the particles 



(see eq.(A.26)) and Lq is the typical angular momentum of 
a generic orbit where the transverse and radial velocities 
are of the same order (see eq.( A.27 )). Besides, the analysis 
detailed in App.JAJ is very general. Indeed, it does not rely 
on the shape of the equilibrium density profile po nor on 
the distribution function /p. As a consequence, this radial- 
orbit instability holds as long as p ^ Lq. This implies that 
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the halo eventually reaches an equilibrium state where the 
transverse velocity v± is of the same order as the radial 
velocity Vr, that is the system becomes roughly isotropic. 
Moreover, this relaxation is very fast since the growth rate 
uj diverges for fj, ^ 0. Indeed, from eq.(p3) the typical 
angular momentum fi of the particles grows as: 



'dt 



ilo M 6 



\/Lo/fi Hot 



(86) 



This merely expresses the fact that the angular momen- 
tum of the particles increases with the time-dependent 
perturbed gravitational potential $i(<). In eq.(|8^) there 
may be an additional power-law prefactor however this is 
irrelevant since the physics is governed by the exponential 
term which is given by eq.dSq). On a small time- interval 
At <^ tu, where to '^ 1/flo is the dynamical time, we 
have t ^to and the growth of the angular momentum is 
well described by: 



^^^oMeVW^ 



(87) 



Then, the time it takes for the angular momentum to grow 
from 0+ up to ALq, with A ^ 1, is: 



T(0^ 



ALn 



1 



dy 
f^o Jo y 



-i/v^ 



This integral converges (very fast) for y ^ 0, hence it 
takes a finite time (T ~ l/^o) in order to go from /i = 0+ 
up to /x ~ Lq. This implies that within one dynamical time 
the typical angular momentum reaches values of order Lqj 
whatever close to exactly radial the system starts from. 
Note that this is very different from the usual power-law 
growth of density fluctuations in the expanding universe, 
where at a finite time to the perturbations can be made 
small enough by starting with a system which is suffi- 
ciently close to uniform. By contrast, from eq.(|88|) we see 
that the system seen after one dynamical time is roughly 
isotropic (f_L ~ Wr), whatever small (but non-zero) the ini- 
tial transverse velocities are. This implies that the func- 
tional pr[5l\ is not continuous at the spherical saddle- 
point derived in Sect. 4.2, Indeed, an isotropic velocity dis- 
tribution provides additional support against the pull from 
the potential well. This means that the halo is somewhat 
more extended than the purely spherical radial solution 
would suggest. This is especially true for power-spectra 
with n < —1 where the radial solution cannot stabilize 
and leads to a slow adiabatic growth of the density as 
particles steadily sinks towards the center of the potential 
well. By contrast, the transverse velocity of the particles 
stabilizes the density profile and the typical radius of each 



particle. For instance, Teyssier et al. (1997) find that in 
the case of spherical gas collapse (i.e. a strongly collision- 
less fluid with an isotropic pressure) the density profile 
stabilizes down to e > 1/6 while for e < 1/6 the isotropic 
pressure is insufficient to stabilize the halo which exhibits 
a density profile of the form p(r,t) ex i(4-24e)/(i8e)^-i^ 
compare with eq. (|79|) and eq.(pO[). We shall come back to 



Therefore, we cannot apply the steepest-descent 
method around the spherical saddle-point obtained in 
Sect. 4.2 to the path-integral dl3). Indeed, as explained 
above, since the functional /^^[(^l] is not continuous at 
this point the spherical dynamics only applies to exactly 
spherical (hence radial) linear density fields which only 
form a subset of vanishing measure. Then, the value of 
the functional pr[5l\ for these spherical states does not 
govern the path- integral (|l3). In simpler words, all real- 
istic density fields show some non-zero deviations from 
exact spherical symmetry which implies, as we proved in 
App.|A|, that the non-linear objects which form in such an 
environment are not described by the known solution of 
the exactly spherical dynamics. 

On the other hand, we note that the analysis detailed 
in App.^ shows that collapsed halos quickly "virialize" , 
in the sense that within a dynamical time their velocity 
distribution becomes roughly isotropic. This is somewhat 
similar to "violent relaxation": starting from an initial 
state which is very far from thermodynamical equilibrium 
(the transverse velocity dispersion (j± is zero) the system 
undergoes a very fast relaxation phase (over one dynami- 
cal time) to reach a new equilibrium state where i7_l is of 
the order of the radial velocity dispersion Ur- This agrees 
with the results of numerical simulations which show that 
within one fifth of the virial radius the velocity field is 



roughly isotropic (e.g., Tormen et al. (1997)). The virial 



radius marks the transition between outer infalling shells 
(which have just experienc ed their first turn-aro und) and 
inner relaxed regions (e.g.. Cole fc Lacey (1996) ). 



this point in Sect. 4.4 



4.4. Virialized halos 

We have shown in the previous section that the minimum 
of the action S[Sl] at the spherical saddle-point obtained 
in Sect. 12 does not govern the path-integral ([l2|). Hence 
it cannot be used to derive the high-density tail of the 
pdf V{pr). Thus, in order to estimate 'P{pr) we need 
to know the behaviour of the functional pr[5l\ for non- 
spherical states. In particular, we are interested in the 
action Sc[Sl\ and the functional pr^c[^l\ defined by con- 
tinuation at spherical points 5l- That is, pr^c[5l\ at such a 
spherical point is defined by the value obtained for a state 
with an infinitesimal deviation A(5l from spherical sym- 
metry in the limit ASl — > 0. Then, we can still expect this 
action S'c[<5l] to show a minimum (or a saddle-point) at 
a spherical state. Indeed, as the deviation from spherical 
symmetry increases the non-linear halo should be more ex- 
tended with an irregular boundary so that a larger fraction 
of the matter lies outside of the radius R. Then, it may 
be justified to apply the steepest-descent method around 
this point. In case one cannot define a continuous action 
5'c[(Jl] (i.e. the hmit pr[5l\ depends on the direction along 
which one approaches exact spherical symmetry) we can 
still expect it can be approximated with a reasonable ac- 
curacy by a smooth functional. This simply means that 
high-density fluctuations probed by V{pr) should be dom- 



p. Valageas: Dynamics of gravitational clustering IV. The probability distribution of rare events. 



15 



inated by such nearly spherical linear density states which 
yield this overdensity pR within the radius R. However, in 
order to be valid this picture requires that these roughly 
isotropic spherical halos are stable (or only weakly un- 
stable). This ensures that a sufficiently large region of 
phase-space (i.e. a subset of initial conditions of non-zero 
measure) is governed by this peculiar dynamics, so that 
it is indeed relevant to the formation of massive collapsed 
halos. 



As we showed in Sect. 4.3, almost spherical non- linear 



density fluctuations become roughly isotropic within a dy- 
namical time. Therefore, since we are interested here in the 
high-density tail of the pdf V{pr<), that is in halos with 
a dynamical time which is smaller than the Hubble time, 
we consider collapsed objects which have already relaxed. 
Thus, we merely need to investigate the s tabi lity of such 
relaxed isotropic halos. As noticed in Sect, 4.3, the trans- 
verse velocity dispersion <j±^ provides some additional sup- 
port against the gravitational attraction hence the orbits 
should stabilize to a finite radius after a few dynamical 
times. Therefore, we expect the stabilized behaviour ( [79| ) 
to apply down to e ^ 0, that is n ^ —3. We shall come 
back to this point below. 

Note that the validity of the density profiles (^9|) rests 
on the assumption that although the velocities quickly re- 
lax to an isotropic distribution the order of magnitude 
of the energy of most particles remains unchanged dur- 
ing the transition phase. We did not obtain a rigorous 
proof of this feature (since the relaxation occurs on a dy- 
namical time the evolution is not adiabatic) but it seems 
quite reasonable. Indeed, the picture we have in mind is 
that inner shells have already relaxed when new shells fall 
in (since the dynamical time associated with each shell 
scales as the time of first turn-around). Then, the new 
particles quickly acquire a significant transverse velocity 
during their first infall so that they "stabilize" onto or- 
bits with a characteristic radius which is of order of their 
turn-around radius. This implies that the density profile 
of the inner regions is not significantly changed and that 
the energy of these outer particles remains roughly con- 
stant. This is to be expected since there is no other en- 
ergy scale in the problem (for a nearly power-law initial 
condition the dynamics is roughly self-similar). This also 
agrees with numerical studies of '"violent relaxation" in 
other contexts ( Kandrup ct al. (1993)[ ). 

Thus, we now consider isotropic halos of radius Re with 
the density profile: 



P{r) = p. ( - 



with 



9e 



3(71 + 3) 



l + 3e 



(89) 



since from eq. (|78|) and eq.(|l9|) we have e = (n + 3)/3. Note 
that for power-spectra of cosmological interest which obey 
—3 < rt < 1 we have < a < 12/5. Of course, the ra- 
dius Re grows with time as new shells turn-around but 
this is not important as the mass within a given radius is 
not governed by the outer shells, as we shall check below. 
Besides, the density profile of the saddle-point is actually 
flat within Ri^ but this is irrelevant here since we only want 



to check whether the overdensity pR is stable with respect 
to the collapse of the outer shells. The gravitational poten- 
tial is obtained from eq. (B9) through the Poisson equation 
which yields: 

*M-*'(i)"° With *.^p4f|5L_ 00) 

for r < Re- Note that with this normalization the value 
of the gravitational potential at infinity is some constant 
*&oo which is not zero, but it is irrelevant for our purposes. 
The shape of <i>(r) is displayed in Fig.pl Thus, we see that 
the cases n < ~1 and n > ~1 are again rather different, 
as could be expected. 



-l<n<l 

2<a<12/5 

Re 




Fig. 2. The shape of the gravitational potential <I>(r) for 
< r < Re. For n > — 1 we have $(r) —^ — oo in the limit 
r — > while for n < — 1 we have <i>(r) — > 0. The plot on the 
right shown for the case < a < 2 actually corresponds 
to 1 < a < 2. For < a < 1 we have the same behaviour 
(i.e <i>(r) -^ for r ^ 0) but the curve looks like r^ rather 
than ^yr. 

Next, we need the distribution function /(r,v) which 
describes these collapsed halos. Indeed, our goal is to show 
that such distributions exist and that they are stable. 
Since the velocity distribution is isotropic we know that 
the distribution function only depends on the energy E of 
the particles (e.g., Binney & Tremaine (1987)): / = f{E). 



Then, the function f{E) can be obtained from the density 
profile p{r) as follows, see Binney & Tremaine (1987). The 
energy E of the particles decreases for smaller r as they 
are more strongly bound to the potential well. Moreover, 
for a halo of size Re the energy of the particles which reach 
the radius Re is E ~ $(i?c) = ^c- Indeed, if they had a 
larger energy they would have a non-zero velocity at this 
point which would imply that some particles can move 
outside of Re, since the velocity distribution is isotropic. 
Then, we define the new variables: 



i/;(r) = $c-$(r) > 0, £ = ^e-E 



^(r)-y> 0.(91) 



Thus, we write the distribution function f{E) as f{£), 
and we have /(£") = for £ < 0. Besides, the density p{r) 
is obtained from /(r, v) by: 



p{r) 



d^v /(r, v) = An 



ip(r) 



d£^2{ij-£)fi£).{92) 
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Then, since both p{r) and tpi''') ^re monotonic functions 
of r we can regard p as a function of ip. Next, differen- 
tiating both sides of eq.(]92) with respect to -0 we ob- 
tain an Abel integral equation which can be inverted as 
Oinney fc Trcmaine (1987)| ): 
£ 



fi£) = 



1 d 



di' 



V^ 



From eq.(B9^ and eq.(pO| 



dp 

d^ 



a 



2-a $, 
which yields: 






1 d£_ 

Ip dip' 

we obtain: 

-2/{2-a) 



(93) 



(94) 



d 



d£ 



2-a $ 

£ 

dv- 



1 



v^^v^ 



1-^ 



-2/(2-a) 



(95) 



We are only interested in the behaviour of f{£) which is 
relevant for small radii r <^ Re where the halo is fully 
relaxed. Indeed, the cutoff of the density profile at the 
boundary R^ is not exact as we should include the tran- 
sition towards the outer shells which are still falling in. 
However, here we only investigate the stability of the in- 
ner regions which are described by power-law behaviours. 
Let us note Is the integral over ^ which appears in the 
r.h.s. in eq.(95). In order to evaluate Is we must sepa- 
rate the cases n > — 1 and n < — 1 which exhibit different 
behaviours. 

First, in the case n > — 1 the energy E of the particles 
which orbit in the inner regions of the halo goes to — oo 
as r — > 0, following the divergence of the gravitational 
potential (<& — > — oo). Therefore, we are interested in the 
behaviour of Is for £ -^ -l-oo, see eq.(pi|). Making the 
change of variable tp — £t we obtain for £ —> +oo; 



Is 



^1/2 



_£_ 



2/(a-2) 



dt 



^2/(a-2) 



(96) 



/o VT-t 
where the integral over t is simply the Euler Beta func- 



tion B(l/2,a/(a-2)), see 
Thus, we get for large £: 



Gradshtcyn fc Ryzhik (1965; 



n > 



with: 



m^N^T:rkH T:ri 



$13/2 



£_ 
1$. 



(6-a)/[2(a-2)] 



(97) 



N, 



1 a{a + 2) 



B 



V8^2 2(a-2)2 
Going back to the variable E we can write eq.(97) as: 
-l<7i<l, 2<a< 12/5, £;< 



(98) 



f{E) ~ TVa 



1$ 13/2 I 1$^ 



^E 



(6-Q)/[2(a- 



■2)] 



(99) 



This gives the behaviour of the distribution function f{E) 
for large negative E, which corresponds to the particles 
orbiting well inside the halo (i.e. r ^ Re). 



Second, in the case n < — 1 the energy E and the 
gravitational potential <I> of the particles which orbit close 
to the center of the halo vanish. Thus, wc arc interested 
in the behaviour of Is in the limit £ -^ $J. Using the 
changes of variables ip = £t and £ = ^cX we obtain for 
x^ 1": 



Is 



$1/2 



dt 



VT^ 



(l-a:t)-2/(2-"). 



(100) 



The integral in the r.h.s. can be written in 
terms of Gauss' Hypergeo metric function 2F1 (see 
iGradshteyn fc Ryzhik (1965)1 §3.197.3). Then, using the 
asymptotic behaviour of the Hypergeometric function for 
a; ^- 1~ we get: 



n<-l: f{£)~M^- 



$ 



3/2 



£_ 



-(6-a)/[2(2-a)] 



with: 



1 0^(0 + 2)^/1 2-fa 



V87r2 2(2 -a)2 \2' 2[2 - a) ^ 
This yields for the distribution function f{E): 



(101) 



(102) 



-3 < n < 



1, < a < 2, < S < $c : 

^ \ -i6-a)/[2{2-a)] 



Pc 



f{E) ^ Af„-9- - 



$ 



3/2 



(103) 



This describes again the orbits which are located within 
the inner regions of the halo. Note that for ?i > — 1 the 
energy spans the range ] — 00, — |<I>c|] while for n < — 1 it 
is restricted to [0, $c]- 

Firstly, we can see that in both cases the distribution 
functions f{E) we obtained are positive since Na > and 
Ma > 0. Hence they correspond to realistic and physi- 
cal distributions. Secondly, we can check that the density 
p{r) at radius r is dominated by the "local" particles with 
an orbit of size '^ r and not by the outer shells of ra- 
dius ^ Re- Thirdly, we can see from eq.(|99|) and eq.(103) 
that in both cases the distribution function / is a decreas- 
ing function of the energy, that is we have d//d£' < 0, 
over the range where / is not zero. Then, we know 
that the property df /dE < ensures that the isotropic 
equilibrium distribution f{E) is stable, see for instance 
[Binney fc Tremaine (1987)| or [Kandrup fc Sygnet (1985)] . 
Therefore, we can conclude that the density profiles (|89[) 
are stable for —3 < n < 1. Thus, after their first turn- 
around the collapsing shells quickly acquire a transverse 
velocity dispersion (7± of the same order as the radial 
velocity dispersion and the particles stabilize to a fi- 
nite radius. This leads to the density profiles (p9) which 
agree with stable-clustering and the collapse of new outer 
shells no longer leads to a slow "compression" of the 
inner mass shells, even in the ra nge —3 < n < — 1. 



Here it is interesting to note that [Teyssier et al. (1997) 
found that the density profiles ( |89| ) only hold for a > 1 
in the case of gas collapse. Therefore, collisionless halos 
are more stable than their hydrodynamical (i.e. gaseous) 
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counterparts. This is actually a rather general fact (e.g., 
Binncy & Tremainc (1987)). Note that the physical pro- 
cesses at work are quite different. Indeed, for a gaseous 
system particles are supported against gravity by their nu- 
merous collisions with neighbouring particles (which gives 
rise to the isotropic pressure) while in a coUisionless halo 
each particle is "stabilized" within the potential well by 
its own kinetic energy. 

4.5. High density tail of the pdfV{pR) 

As explained in the previous sections, the high-density tail 
of the pdf V{pb) should be governed by almost spherical 
halos which have relaxed to a roughly isotropic equilib- 
rium distribution function. The density profile is stabi- 
lized by the transverse velocity dispersion so that stable 
clustering holds. Moreover, these equilibrium states are 
stable (at linear order) with respect to small perturba- 
tions. Therefore, we can now apply to the whole range 
—3 < n < 1 the considerations used in Sect. 4.2 for ti > —1. 



More precisely, we approximate the functional pr[6l.r"] 
by eq.(17) where the Lagrangian scale Rl is associated 



with the particles with a turn-around radius equal to R. 
The function !Fp{5l) now describes the continuous limit of 
almost exactly spherical collapse. In other words, it is not 
given by the exact spherical dynamics with radial trajecto- 
ries: it is defined as the non-linear overdensity pn obtained 
in the limit of infinitesimal deviations from spherical sym- 
metry. Thus, it corresponds to the virialized distribution 



functions f{E) obtained in Sect. 4.4 



As recalled in Sect. 3. 1 , a relation of the form (h/n) leads 
to a spherical saddle-point given by eq.(p^ and eq.(p^. 
Provided the steepest-descent approximation is justified 
in the limit we consider here (i.e. pn — > oo) this yields for 
the generating function "0(2/) the expression (p3), where 
eq.(24) through eq.(^) must be used with the relevant 



function J-p{5i,). We refer again the reader to paper II for 



Tp{5l) = (1 + a,) ( ^ 



(104) 



Note that this is very similar to the usual Press-Schechter 
prescription (Press & Schechter (1974)). In particular, the 
normalization parameters Ac and 5c should be close to the 
usual values 5c — 1-69 and Ac ~ 177. However, these lat- 
ter values are only approximate. In order to obtain the 
right normalization one should run a numerical simula- 
tion in order to study the collapse of a density profile 
of the form ([l5[), to which is added a small deviation 
from spherical symmetry. In other words, one must ob- 
tain the behaviour of the saddle-point of the continuous 
action S'c[(5l]. However, such a numerical study is beyond 



the scope of this paper. Note that we can expect a small 
dependence on the slope n of the power-spectrum of the 
normalization (1 + Ac)/(5c. From eq.(p9|) we obtain the 
characteristic function Qp{t) as: 

2/(n+5) 



r«-l: gp{T) 



1 + Ac 
5-^ 



i-r) 



6/(n+5) 



(105) 



Note tha t large o verdensities p^ are associated with t -^ 
— oo, see paper II . Then, from cq.(|25|) we have the asymp- 
totic behaviour for pn. — > +oo: 

/1 + A \2/("+2) / -6 z/ \ ^/^"+'^ 



Note that for n > —2 large positive pn are associated with 
large negative y. This corresponds to a pdf V{pr) which 
decreases faster than a pure exponential at large densities. 
By contrast, for n < — 2 high densities are associated with 
y — > 0^ , which gives rise to a cutoff which is smoother than 
an exponential. This is discussed at length in [paper II , 
where we come across the same transition at n = in the 
quasi-linear regime. 

As for the case of underdensities described in Sect. 3.2 
we now need to check whether the steepest-descent 
method is justified in the limit of high densities. This is 
again a difficult point since we do not know the exact form 
of the functional pu[5l\- However, we can apply again the 
discussion developed in Sect. 3.3 about the ordinary in- 
tegral (p3) since we have similar power-law behaviours. 
Then, the generating function tp{y) is still given by the ex- 
pression ( p6| ) and as in Sect. 3.4 we use the approximation 
(BTI) from the analogy with the one-dimensional case, tak- 
ing into account the expansion of the universe. Note that 
the factor l/p_R is now associated with "sharp" directions 
(the translations). Then, we obtain again the expression 
( 461) for the pdf V{pfj). Using the power-law behaviour 
( |l05|) this yields: 

1/3 



a detailed derivation of these points. Thus, we now need to V{pr) 
write an explicit expression for the function J-'p{5l)- Since 
after virialization we assume the density profile to remain 
stable the overdensity pn grows as p^ ex p~^ ex a"^ while 
5l{M) cx a (here we consider a critical density universe). 
Hence we get pR oc 5\ and we write: 



1 



27rcr 6 



1 + Ac 
<53 



Pr 



(7-n)/6 



-(^ 



,)-2/3^(. + 5)/3/(2^2) 



(107) 



As noticed above, for n < — 2 the saddle-point is actually 
a local maximum of the action so that the steepest-descent 
method requires some care. This is similar to what occurs 
for n < in the quasi-linear regime. Thus we refer the 
reader to App.A in paper 11 for a detailed discussion of 



the way to handle such cases. However, eq.( |l07| ) remains 
valid for all n. 

Of course, as for the case of underdensities discussed 



in Sect. 3.6.2 the expression (107) agrees with the sim- 
ple spherical model (p^), where the spherical dynamics 
used to relate the linear and non-linear density contrasts 
is given by eq.(|l04|). Thus, we rec over the pdf obt ained 
in Valageas fc Schaeffcr (1997) and Valageas (1998) from 
the spherical model couple d w ith the stable-clustering 
ansatz. As di scus sed in Sect, 3^, the prefactor which ap- 
pears in eq.( |l07| ) is not exact because we used a sim- 
ple approximation for the determinant T) which arises 
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from the Gaussian integration of the path-integral around 
the spherical saddle-point. In particular, we can expect 
the exact calculation to give a numerical multiplica- 
tive factor of order unity. However, we note that the 
same approach used in the quasi-linear regime provides 
very good results as compared to numerical simulations 
(e.g., Valageas (1998), paper II). Therefore, the correction 
might be quite small (but this is a different regime hence 
we cannot draw definite conclusions). 

The pdiVipR) shows a cutoff at the characteristic den- 
sity pc such that Tc ^ —a. In the quasi-linear regime this 
corresponds to Sr ^ a while in the non-linear regime 



where we can use eq.(105) this yields p^ ~ a 



T-6/(n+5) 



Therefore, we can write the cutoff pc as: 



Pc = Min l + cr(i?),cr(i?) 



6/(n+5) 



(108) 



which is valid for all values of a. Indeed, we recall that 
we only used the limit of rare overdensities and the value 
of (T is actually irrelevant, so that our results apply from 
the quasi- linear regime (a <^ 1) up to the highly non- 
linear regime {a ^ 1). The influence of cr(-R) only appears 
in the cutoff pc. As we probe deeper into the non- linear 
regime the characteristic overdensity pc of virialized halos 
is pushed to higher values. In the highly non-linear regime, 
eq.(108) gives: 



cr(i?)>l: Pa = (TiR)~ 



3(>i + 3) 
,1 + 5 



aR) (109) 



where we noted Ro the non-linear scale (defined by 
cr(i?o) — 1) and we used eq.(p8|) for the last term. Here, 



we must point out that we only derived eq.( |l07| ) in the 
limit of rare overdensities. That is, it only applies to large 
densities with pn ^ pc- The behaviour of the pdf 'P{pr) 
at lower densities cannot be obtained by the steepest- 
descent method detailed in this article. Thus, for inter- 
mediate densities pv < pR < pc the path-integral (|lj) is 
no longer dominated by one saddle-point. Then, one must 
devise another approach which takes into account the dy- 
namics of a wide range of initial conditions Sl (x) . Going 
back to the interpretation of the results (E^) and (107) 



in terms of the simple spherical model this is also quite 
clear. Indeed, these intermediate densities correspond to 
relatively low densities (i.e. smaller than the typical den- 
sity Pc ~ C associated with collapsed halos of size R) where 
virialization processes occured lately. Then, the dynamics 
of these objects must have been strongly influenced by 
the gravitational interaction (e.g., mergings, tidal effects) 
with neighbouring halos which are typically more massive. 

4. 6. Aging processes 7 

The very dense halos with pji ^ pc are sufficiently 
rare and massive not to be affected by the interaction 
with neighbours as they form. However, we note that the 
present context shows some important differences with the 
case of rare underdensities studied in Sect.||. Indeed, in 



that former case the underdensities which govern the low- 
density tail of the pdf consist of rare regions which are still 
expanding in the universe and the dynamics is fully de- 
termined by the spherical model. Moreover, outer regions 
have no strong influence on the behaviour of the under- 
density. By contrast, in the present case virialized halos 
remain stable after their collapse. However, this may only 
be a zeroth order approximation. Indeed, as larger scales 
turn non-linear and collapse the halo of size R which we 
study becomes embedded within more massive and ex- 
tended structures. Then, as time goes on we can expect 
repeated interactions with neighbouring halos and tidal 
effects (e.g., mergers, variations of the large scale gravi- 
tational potential) to have a cumulative impact onto this 
object which will gradually be distorted or even disrupted 
or merged within a larger halo. Note that these "aging 
processes" should be less efficient for more massive ob- 
jects. Therefore, a possible scenari o wo uld be that the 
high-density tail we obtained in eq.(107) only applies to 
densities which are larger than a time-dependent thresh- 
old pt which increases faster than p^ with time. Thus, we 
would have pt ~ Pc for ^ ^ 1, when larger scales have 
not collapsed yet, and pt/pc —^ oo for ^ — *■ cxi, when tidal 
effects (in a broad sense) have had plenty of time to influ- 
ence the density profile of the halo. 

Let us briefly describe how such a mechanism might 
show up in the path-integral formalism we used in the 
previous sections. As seen from eq.dl2) and eq.(O), the 
term (J^.A^ .6l) in the exponent of the path-integral 
leads to a characteristic cutoff 1(5^.^/1 ^ cr(i?') for the non- 
spherical perturbations at scale R' of the density fields 5l 
around the saddle-point which provide a significant con- 
tribution to tpiy)- These deviations are much smaller than 
the density Sl.b/ of the saddle-point, in the limit of rare 
events, for R' ^ R^. Therefore, they have no strong im- 
pact on the collapse of the halo. However, on much larger 
scales they become non-negligible since the density profile 
of the saddle-point declines much faster at large scales, 
as (Sl.a' oc a{R'Y as seen in eq.([l9|). Then, when these 
scales turn non-linear (i.e. much after the halo of Eulerian 
size R has collapsed and virialized) these important de- 
viations from spherical symmetry may play a key role. 
Indeed, since the collapse is no longer roughly self-similar 
but proceeds in a rather irregular manner, because of the 
break of spherical symmetry, one expects the formation of 
an irregular distribution of clumps. Then, this could lead 
to the tidal effects or merging processes described above. 
Indeed, the regular and adiabatic collapse of outer shells 
(adiabatic with respect to the inner shells) is replaced by 
an irregular evolution which may proceed through sudden 
non-linear changes in the large scale gravitational poten- 
tial. This means that the region of phase-space where the 
action S[5l] is close to its value at the saddle-point be- 
comes narrower as time goes on: one needs to restrict to 
linear states 5l which are increasingly spherically sym- 
metric. This implies a change in the normalization of the 
generating function ^(y) in this regime. Mathematically, 
this arises from an increase of the determinant T) which ap- 



p. Valageas: Dynamics of gravitational clustering IV. The probability distribution of rare events. 



19 



pears in eq. (pq) (transposed to the regime of high overden- 
sities) from the Gaussian integration around the saddle- 
point. In other words, the action S[Sl\ exhibits an increas- 
ingly fast variation with Sl as time goes on. 

The effect of these processes on the high-density tail 
of the pdf V{pr) is now clear. The exponential cutoff ob- 
tained in eq.(107) remains valid (at least in a first step) 
but the prefactor steadily decreases as time goes on, fol- 
lowing the decline of the phase-space volume (i.e. the num- 
ber of states 5l) which leads to this non- linear overden- 
sity PR. Note however that this decline could be hidden if 
these tidal effects also build deeper minima of the action. 
These effects should be stronger for lower overdensities pR. 
Finally, after the pdf 7'(p/i;) at such a point has undergone 
a significant evolution we can expect that it enters a new 
regime (possibly non-stationary) where it is governed by 
these "tidal" processes. It is clear that this regime cannot 
be obtained by a steepest-descent method similar to the 
approach developed in this paper. We must also note that 
these arguments, although quite reasonable, are still hy- 
pothetical at this stage for want of a rigorous theoretical 
derivation. However, this point requires new theoretical 
tools and a detailed analysis which are beyond the scope 
of this article. 

4. 7. Comparison with the non-linear scaling model 



Here we briefly compare our resul t (107 ) with th e no n- 
linear scaling model recalled in Sect. 3.6.3 . Since eq.(107) is 
consistent with the stable-clustering ansatz -as it neglects 
the possible "aging processes" discussed in Sect. 4.6- 



it IS 



clear that it is also consistent with the scaling laws (|57|). 
Let us recall briefly the main properties of the non-linear 
scaling model (p7|) with respect to the high-density tail 



of the pdf V{pr). As shown in Balian & Schaeffcr (1989) 
the scaling laws (|5^) imply that for large overdensities 
PR ^ Pv the pdf exhibits the scaling: 



PR > Pv ■ ViPf 






with 



PR 



where we defined the scaling function h{x) by: 
/■+*°° dw 



(110) 



(111) 



and the function ip{y) was introduced in eq.(p2f). As shown 
Schaeffcr (1997)[ the expression ( |107 



in Valageas 

actually be recast in the form ( |110D . Indeed, for a power 
law linear power-spectrum P{k) ex k"' the non-linear cor 
relation ^{R,t) can be written as a function of a{R,t). 
From cq.dSq) we obtain S, c>c (t6/("+^). Therefore, we write: 



a > 1 : S,{R) 



1 



<53 



2/(n+5) 



,(i^)6/(«H 



5) 



where both correlations ^ and a are taken at the same 
Eul erian scale R. We wrote the normalization factor in 
eq.(112) as in eq.(105) and 1 -I- A^ is an unknown param- 
eter which may be taken from numerical simulations. It 



corresponds to an "effective" non-linear density contrast 
associated with the linear contrast 5c. This agrees with 



the scaling ansatz introduced by Hamilton et al. (1991) 
to relate the non-linear correlation at scale R to the linear 
variance a{RL) at the Lagrangian scale R^. In particular, 
as noticed in Valageas & Schaeffcr (1997) (App.E) the re- 



lation (j{Rl) '^ S.{R) is well described by the spherical col- 
lapse model. Thus, eq.(112) actually means that we have 
£,{R) = J-^[a{RL)] where JF^ is of the form (104) where 
Ac is replaced by A^. Note that we should have A^ < A^ 
since density fields with important deviations from spher- 
ical symmetry should reach a smaller non-linear density 
contrast A than the value Ac obtained for the (almost) 
spherical saddle-point and A^ corresponds to a n av erage 
ove r all possible density fields. Sub stitu ting eq. (|ll2| ) into 
eq.(107) we obtain the scaling law (110) with: 



hs{x) = 



1 



/2n 6A 



a.-(7-n)/6 g-:E("+^'/V(2A^) 



(113) 



where the subscript "s" refers to the spherical model, and 
we defined the parameter A by: 



A = 



1 



1 



1/3 



^1. 



(114) 



The parameter A is expected to show some dependence on 
n since both Ac and A^ depend on n. This dependence 
for Aj was indeed checked in numerical simulations, see 



for instance Jain et al. (1995). 

As we pointed out in Sect.[1.5| the expression (107) only 
applies to rare overdensities pr^ Pc- In terms of the scal- 
ing variable x introduced in eq.(|llC| ) this means that the 
function hg{x) written in eq.( |113| ) only applies to x 3> 1. 
Unfortunately, this makes the comparison with numerical 
simulations rather difficult since one needs to probe the far 
high-density tail of the pdf 7'(/9jj). Indeed, the points mea- 
sured in current numerical simulations from counts-in-cells 
statistics do not g o much beyond x ^ 10, see for instance 
Fig. 4 and Fig. 6 in Valageas et al. (2000) This is not suf- 
ficient to really probe the cutoff of the pdf. However, in 
Valageas et al. (2000) we devised a statistics which allows 
one to go slightly deeper into the high-density tail. Thus, 
instead of computing the pdf V{pii) one can investigate 
the cumulative mass function Fr{> Af ) of halos of radius 
R more massive than M. In practice, one uses a "spherical 
overdensity algorithm" , looking at particles in order of de- 
creasing density and defining halos as objects of constant 
size R around the density peaks. Then, the differential 
mass function iiR{M)dM/M should obey the scaling law: 

dAf T„, , da; 



MM)- 



M 



'■H{x) 



(115) 



with a scaling function H{x) which is close to h{x), see 



(112) [Valageas fc Schaeffcr (1997)| and [Valageas et al. (2000) 



In particular, the high-density tails are expected to show 
the same beha viour (up to a normalization factor of or- 
der unity), see Valageas fc Schaeffcr (1997) . This proce- 



dure allows one to probe slightly deeper into the high- 
density tail of h{x) because "cells" are directly drawn 



20 



P. Valageas: Dynamics of gravitational clustering IV. The probability distribution of rare events. 





' ' ' ' ' 











_ 










n-O 


- 






. 


^'^r^. 


- 


- 




^^^^ 


f^^^ 


P§^ 


k 


:^ 


j:.^^ 


W^ 


/ 










□ -1^ - 


0,5 




\ 


- 


X 


X R = 
_ R = 


1 






" 


y^ 


/x R = 


4 




\ ; 


- 




1 , , 


, , 1 




1 



that the numerical simulations arc consistent with the ex- 



log(x) 



_ 







1 




_ 


- 






n-- 


- 1 


- 


- 


^ 








- 


- 






- 


- 


x"^ 




''*A 




- 




/^ 




. ^A 






_ 


X^ □ 1-; - o,o 




. -^ 




_ 


- 


/^ X R - 1 
/^ _ R = 2 




J 




- 


- 


^ ^R = 4 






\ 


- 


-y^ 










- 


- 


1 , , , , 1 , , , , 




1 


1 


- 




log(x) 



Fig. 3. The mass functions of hales defined by various co- 
moving radii i?, obtained from a "modified spherical over- 
density" algorithm. The data points are taken from the 



numerical simulations described in Valageas et al. (2000) 



Different symbols correspond to different values of R. Note 
that for each radius we display the results obtained at 
several times which correspond to various values of ^. The 
solid curve is the scaling function hs{x). We use the values 
A = 3.2, 2.8 and 2.1 for n = 0, -1 and -2. 



poncntial tail written in eq.(113). Note however that we 
expect the normalization of eq.(113) to be only approxi- 
mate because of the approximation ( pT| ) for the determi- 
nant v. Moreover, as noticed above H{x) may differ from 
h{x) by a factor of order unity, in-between 7/3 and 1 in 
simple models ( Valageas fc Schaeffer (1997)1 ). 

Note that the data points shown in Fig.y are obtained 
from different comoving radii as well as from the same 
comoving radius at different times (hence different val- 
ues of £, and a). Therefore, the fact that all curves su- 
perpose shows that the stable-clustering ansatz ( p7{ ) is a 
good approximation in the regime probed by these nu- 
merical points. In any case, the point of FigJ^ is merely 
to show that eq.(107) is consistent with numerical simu- 
lations. In fact, it is difficult to see how one could "beat" 
the exponential cutoff given by eq. ( |107| ) . Indeed, it is clear 
that in the rare-event limit the cutoff of the pdf is di- 
rectly linked to the initial Gaussian cutoff. Then, in order 
to change the high-density tail obtained in eq.(|l07|) one 



should rely on a mechanism which would violate eq. ([KM 
As explained in Sect. [4.6| such a process may indeed ex- 
ist (through gravitational interactions with outer shells) 
but it should be negligible if one considers sufficiently 
large overdensities pR (at a given time). Of course, we 
can see in Fig.0 that the pdf (107) fails at small densi- 
ties (i.e. X ^ 1). As explained in Sect.E^ this is quite 
natural since the steepest-descent method should not ap- 
ply to this regime. Therefore, in order to derive a predic- 
tion for this low-density part of the pdf (which seems to 
exhibit a power-law behaviour) one needs to build other 
non-perturbative tools which fully take into account the 
involved processes of mergings and tidal interactions (with 
similar or more massive neighbours) which govern the 
dynamics of these intermediate objects. Here we must 
not e th at in Valageas et al. (2000) we had concluded that 
eq.(113) does not agree with the counts- in-cells 7'(A^) mea- 
sured in the simulations. However, we used different values 



of A and we compared the pdf (107) with numerical points 
over the whole range pr » pv Here we still agree with 
this conclusion (as noticed above the small— cc tail clearly 
fails) but we note that the high-density cutoff (a; 2> 1) 
may be consistent with numerical simulations. 



4.8. Press-Schechter mass function 



around high-density peaks. Then, these high-density fluc- 
tuations are well accounted for, while for standard counts- 
in-cells statistics (i.e. the pdf VipR)) such high-density 
peaks are usually split over several cells. Therefore, we 
compare in Fig. ra the scaling function hs{x) obtained in 
eq.(113 ) with the results o f numerical simulations taken 
from Valageas et al. (2000) (note that this corresponds to 
the PS prescription without the factor 2). We choose the 
parameter A so as to get a reasonable fit to the numerical 
points (see the caption for the relevant values). We can 
check that we get A ^ 1 as it should. The figures show 



Finally, in view of the practical importance of the Press- 
Schechter (PS) prescription (Press & Schcchter (1974)) to 
compute the mass function of just-virialized objects we 
comment this "recipe" in the light of the present work. 
The PS approach provides a simple estimate of the cu- 
mulative mass function F{> M) which gives the fraction 
of matter embedded within just-virialized halos (i.e. with 
a non-linear density contrast Ac) of mass larger than M . 
As for the spherical model ( |5l| ) it merely approximates 
F(> M) by the fraction of matter which is enclosed within 
spherical cells of Lagrangian radius Rl and linear density 
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contrast greater than Sl,Rl, with Sl,Rl — ^c- Therefore, 
as in eq.llsl^) it writes: 



16) 



/•OO /'OO 

Fps{> M) = / dSLVLiSL) = / dv V^l\v) (1 



where the subscript "PS" refers to the Press-Schechter 
prescription. This gives the differential mass function 
fiPs{M)dM/M as: 



^lPs{M) 



dM 



dF, 



PS 



/2^ cr(M) 



dM 
diner 



dM 



dlnM 



-Sl/{2aiMy) 



dM 



(117) 



Then, in order to obtain a mass function which is cor- 
rectly normalized to unity it is customary to multiply the 



expression (117) by an ad- hoc factor 2. 

The quantity we investigated in this article is not the 
mass function of just-virialized halos but the pdf Vipp), 
that is the statistics of the counts-in-cells. In order to com- 
pare our results with the PS mass function we first need a 
prescription to obtain the mass function from the counts- 
in-cells. To do so, we can employ the prescription used in 
the PS approach, which we now apply to the non-linear 
density field. That is, we simply write: 



F(> M) 



dSil + 6)Vid). 



(118) 



This means that F{> M) is approximately the fraction 
of matter which is enclosed within spherical cells of ra- 
dius R and linear density contrast greater than Ac. This 
is the translation of eq.(116) in terms of the non- linear 
density field. Note however that in both cases this is only 
an approximation since it does not t ake i nto account the 
cloud- in-cloud problem. The relation ( 118 ) is more g enera l 
and should yield better results since, contrary to eq. (|ll6| ), 
it does not assume that there is a one-to-one relation be- 
tween the linear and non-linear density contrasts (given by 
the spherical dyn amics for the PS prescription) . In partic- 
ular, as shown in Valageas fc Schaeffer (1997) , the cloud- 
in-cloud problem is very severe at small masses when 



one works with the linear density field as in eq.(116) so 
that the low-mass tail of the mass function is actually ill- 
defined. By contrast, simple models (based on the stable- 
clustering ansatz recalled in Sect. 3. 6. 3| ) show that when 
one works with the actual non-linear density field as in 
eq.(118) the cloud-in-cloud problem becomes much less se- 
rious since it only yields a correction of order unity for the 
normalization of the low-mass tail of the mass function. 



see 



Valageas & Schaeffer (1997). Of course, the problem 
with eq.( |118| ) is that one needs to estimate the non-linear 
pdf V{5ii), which is still an unsolved problem (in this ar- 
ticle we only obtained the tails of this p df). 

As noticed in Sect. IS , the expression ( |l07|) agree s with 
the simple spherical mo del (|5l| ) detailed in Sect. 3.6.2 . 



Therefore, using eq.( |ll8D we recover the usual PS pre- 
scription without the factor 2. Her e, w e must note that 
the prefactor which appears in cq.( |l07|) is not exact, as 



discussed above, so that a rigorous calculation of the de- 
terminant V might yield a factor 2. However, this is rather 
unlikely and until wc have a rigorous estimate of this de- 
terminant we can as well keep the normalization obtained 
in cq.(107) or simply fit the normalization to the results of 
numerical simulations. Note that since the PS mass func- 
tion deals with "just-virialized" halos we do not encounter 
the "aging processes" discussed in Sect. 4.6. Indeed, these 



halos have just collapsed and they have not had time to 
suffer from the cumulative effects brought by the interac- 
tion with smaller neighbours. Therefore, the "exponential" 
form of the cutoff of the PS mass function should be exact. 



Wc must stress that the expression (|107[ ) was only de- 
rived in the limit of high overdensities, that is for p/j 3> Pc. 
Thus, as discussed in the previous section, our analysis 
shows that the exponential cutoff of the PS mass function 
is correct but the PS approach should not be used for the 
low-mass (or low-density) tail. In other words, the power- 
law behaviour which appears in eq.( |l07| ) for pp ^ pc, 
or for small masses M ^ M^ in the PS mass function, 
is not justified by the approach developed in this arti- 
cle. In fact, we can actually expect this power-law to be 
wrong since we know that deviations from spherical sym- 
metry play a key role in this regime. Besides, as we recalled 
above the cloud-in-cloud problem removes any predictive 
power to the low-mass tail of the PS mass function (e.g., 
[Valageas fc Schaeffer (1997) ). This implies that the high- 
mass cutoff of the PS mass function, or the high-density 
cutoff of the pdf 'P{pr) written in eq.(107), should not 



necessarily be multiplied by a factor 2 in order to obtain a 
correct normalization of the mass function to unity, since 
the usual analytical mass function (or pdf) does not ex- 
tend to M < Ale or PR < Pc- 

Moreover, the values Sc — 1.69 and Ac = 177 used 
in the literature are only approximate. We did not derive 
the exact values in this paper but we explained in Sect. 4.4 
how they could be obtained by numerical means. One 
needs to find the saddle-point 5l(x) of the action 5[(5l]. A 
simple procedure which should be sufficient for practical 
purposes is to study through a numerical simulation the 
dynamics of a spherical linear state of the form (|l5|), to 
which we add a small perturbation from spherical symme- 
try. Then, one could measure the function Tp{5l) written 
in eq.(104). Of course, we expect to recover values close 
to the standard density contrasts 5c = 1.69 and Ac = 177 
but there should be a small dependence on the slope n of 
the power-spectrum. Here, we note that recent numerical 
simulations (e.g., Governato et al. (1999)| ) show that the 
usual PS mass function (with the factor 2) overestimates 
the number of small halos and underestimates the number 
of massive halos. In the light of the work presented in this 
article, the deviation at small masses (below the cutoff 
Mc) is not surprising since we do not expect the PS mass 
function to apply to this regime. On the other hand, the 
exponential tail of the PS mass function should be correct. 
Then, the discrepancy with the numerical results would be 
due to the fact that these authors (following the standard 
practice) do not use the exact values of the density con- 
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trasts {Sc, Ac) nor the correct normalization. Indeed, as we 
explained above there is no reason to multiply the high- 
mass tail by a factor 2. Unfortunately, we did not manage 
to predict accurately this normalization factor. Moreover, 
there still exists the possibility that an exact calculation 
of the determinant V which appears through the Gaussian 
integration around the saddle-point modifies the exponent 
of the power-law prefactor in eq.(107) or gives rise to log- 
arithmic corrections. 

5. Conclusion 

Using a no n-perturb ative method developed in a previ- 
ous work ( paper U ) we have investigated in this arti- 
cle the tails of the pdf V{pr). Our approach is based 
on a steepest-descent approximation which should yield 
asymptotically exact results in the limit of rare events. 
Therefore, it applies to all values of the rms linear den- 
sity fluctuation a, from the quasi-linear up to the highly 
non-linear regime. 

First, we studied the low-density tail of the pdf, that 
is very rare underdensities. This regime is rather simple 
since shell-crossing only occurs for power-spectra with a 
slope n > — 1, and even in this case shell-crossing only has 
a limited impact as it merely introduces numerical factors 
of order unity in the pdf but it should not modify the char- 
acteristic exponents which govern the low-density cutoff. 
This allows us to get a good description of rare underden- 
sities and to derive the shape of the low-density tail of the 
pdf. Apart for the power-law prefactor (which we did not 
derive in a rigorous manner) the exponential cutoff should 
be exact, provided there are no deeper non-spherical min- 
ima of the action, which appears rather unlikely. Besides, 
we have shown that this regime could be recovered from a 
simple spherical model. This is due to the spherical sym- 
metry of the problem we investigate, which implies that in 
the limit of rare events the pdf is governed by rare almost 
spherical density fluctuations. 

Moreover, we have shown that our results agree with 
perturbative calculations over the range —3 < n < — 1 
where the latter yield finite results (up to the first sublead- 
ing term). Besides, since our method is non-perturbative 
it also applies to — 1 < n < 1 where shell-crossing comes 
into play (which leads to the break-up of perturbative ap- 
proaches). This makes the derivation of the exact low- 
density tail more difficult but we obtained an approximate 
result which should provide the exact exponents which 
govern this low-density falloff. Note that in this regime 
our analysis shows that the spherical model only yields an 
approximate result for the low-density tail. Moreover, it 
clearly shows that perturbative results should be viewed 
with caution since leading-order terms may turn to be 
wrong (even when they are finite). This can be understood 
from the fact that perturbative expansions diverge (in 
fact, beyo nd a finit e order one encounters divergent quan- 
tities, see paper V ). Next, we have pointed out that our 
results show that this low-density tail cannot be obtained 
from the stable-clustering ansatz which is often used to de- 



scribe the highly non-linear regime {a 3> 1). However, this 
does not imply that the latter model is wrong. It simply 
means that it is a zeroth-order approximation which does 
not capture the behaviour of these rare expanding voids. 
Finally, we have noticed that the low-density tail we de- 
scribed in this work is still out of reach of current numer- 
ical simulations. Therefore, because of the finite number 
of particles available in these numerical works, they can- 
not probe the actual low-density cutoff of the underlying 
continuous matter distribution yet. 

Second, we have turned to the high density tail of the 
pdf. This case is much more difficult since shell-crossing 
now plays a key role. We have shown that a naive ap- 
proach based on the exact spherical dynamics (which im- 
plies radial trajectories) fails. Indeed, a strong radial-orbit 
instability implies that the radial collapse solution is ac- 
tually irrelevant. In particular, we have shown that col- 
lapsed halos see their velocity distribution become roughly 
isotropic over one dynamical time, whatever small the ini- 
tial deviations from spherical symmetry. This leads to a 
very efficient virialization process, similar to "violent re- 
laxation" . Starting from an initial state which may be very 
far from thermodynamical equilibrium (in the sense that 
the transverse velocity dispersion may be very small) the 
halo relaxes over one dynamical time to a new equilibrium 
state where the velocity distribution is roughly isotropic. 
Moreover, we have shown that this process stabilizes the 
density profile. Thus, contrary to the cases of radial col- 
lapse or gaseous dynamics, the transverse velocity disper- 
sion stabilizes new infalling shells to a finite radius so that 
the almost spherical halo obeys stable-clustering. 

Then, using these results we have derived the high- 
density tail of the pdf V{pr). The exponential cutoff 
can again be recovered from a simple spherical model 
(apart for a possible modification of the power-law prefac- 
tor). Moreover, this is consistent with the stable-clustering 
ansatz in the highly non-linear regime. We have also shown 
that these results are consistent with numerical simula- 
tions. Besides, it implies that the exponential cutoff of the 
standard Press-Schechter (PS) mass function is correct, 
although the prefactor and the characteristic density con- 
trast 5c may be modified. However, there is no reason to 
multiply the PS prescription by an ad-hoc factor 2 since 
it should not apply to low-mass halos. Similarly, our re- 
sults only give the very high-density tail of the pdf (i.e. 

PR > D- 

In-between these extreme low-density and high-density 
tails of the pdf V{pr) there appears in the non-linear 
regime (cr <; 1) an intermediate region which cannot be 
described by our steepest-descent approach. This range 
disappears in the quasi-linear regime (tr — > 0) as any fi- 
nite density contrast becomes a rare event in this limit. 
This is why our approach can full y descri be the qua si- 
linear regime as detailed in paper II (see also paper III for 
non-Gaussian initial conditions), in a way which is fully 
consistent with the study developed here. By contrast, in 
the non-linear regime this intermediate range corresponds 
to "moderate" density fluctuations whose dynamics is 
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strongly influenced by their neighbours (e.g., tidal effects, 
mergings). Moreover, we argued that this regime might ex- 
tend towards larger densities (in physical coordinates) as 
time goes on and the sensitivity to deviations from spher- 
ical symmetry grows. This would imply that the range of 
validity of the high-density tail we derived in this article 
would be repelled to increasingly large densities as one 
probes deeper into the highly non-linear regime. An un- 
derstanding of this regime requires new non-perturbative 
tools which can handle the intricate non-local dynamics 
of these typical regions through strong tidal effects and 
merging processes. It is not clear at present whether this 
regime can be described by the stable-clustering ansatz, 
although numerical simulations suggest it should provide 
at least a good zeroth-order approximation. 

Finally, we note that although we considered the case 
of a critical-density universe in this paper, it is straight- 
forward to extend our results to arbitrary cosmological 
parameters (rim,riA)- One simply needs to use the rel- 
evant function J^p{5l) which describes the spherical dy- 
namics (taking into account virialization) . This does not 
introduce any new feature. Moreover, although we consid- 
ered power-law linear power-spectra P{k) ex /c" our results 
can easily be extended to other smooth power-spectra like 
the usual CDM model. This merely modifies the char- 
acteristic function Qpir) defined from the spherical dy- 
namics function J-p{Sl), through the factor a{RL)/o'{R) 
which is no longer a power of the overdensity. However, 
for practical purposes it should be sufficient to use the ex- 
pressions derived here for power-law power-spectra, sub- 
stituting for the index n the local slope of P{k) over 
the range o f wavenumbers one is stu dying, as done for 
instance in Peacock fc Dodds (1996) to obtain the non- 



linear two-point correlation function from its linear coun- 
terpart. Besides, it is clear that our results can also be ex- 
tended to non -Gaussian initial conditions, like the models 
we studied in paper III . In all cases, the simplest way to 
obtain the relevant expressions is to use the simple spher- 
ical model recalled in this article coupled with the pdf of 
the linear density field, in a fashion similar to the Press- 
Schechter prescription. This shows no difficulties. The only 
technical point is to derive the relevant linear pdf, but this 
can be done in a rigorous manner following the results de- 
tailed in paper III , or through an appropriate application 
of the steepest-descent approach developed in paper 11 to 
specific models. 

Appendix A: Radial-orbit instability 

In this appendix, we investigate the stability of a spheri- 
cally symmetric halo where all particles follow radial tra- 
jectories. Such a system arises from the collapse of an ini- 
tial density fiuctuation which obeys exact spherical sym- 
metry for the linear growing mode. As expected, we shall 
find below that such a system exhibits a very strong radial- 
orbit instability with a growing rate u> which goes to in- 
finity as the trajectories become closer to radial. Thus, for 
our purposes we can restrict ourselves to small time in- 



tervals At which are much smaller than the Hubble time 
tn- Then, we can neglect the growth of the halo as new 
outer shells collapse on the time-scale tn- Therefore, we 
study here the stability of a stationary spherical halo. If 
all orbits were exactly radial the equilibrium distribution 
function /o(r, v) would be of the form g{E)5o{L'^), where 
E is the energy of the particle and L = |L| the magnitude 
of the angular momentum. However, since perturbation 
theory is singular around such a state (as we shall see be- 
low) we add a small angular momentum to the particles. 
Thus, we write: 



/o(r,v) 
with: 



g{E)K{L) 



(A.l) 



/ di^ hp{L) = I dL^ Sd{L^) = 1 (A.2) 

and hp{L) is peaked for very small values of L of order fi. 
For instance, we write: 



^^(^) - ^ K^ 



with 



da;^ h{x) = 1, 



(A.3) 



and we investigate the limit /i ^ 0+ . We use the spherical 
coordinates (r, 9, 0; Wr, f j_, a) with: 



vj_ = Jvg + v-^, vg = v±_ cos a, W0 — v±_ sin a. (A. 4) 

The angle a spans the range < a < 2%. We note $o(^) 
the equilibrium gravitational potential and (/i,$i) the 
perturbed distribution and potential at linear order. The 
linearized collisionless Boltzmann equation reads: 

dfo 



— +v.V-V$o.T^)/i = V$i. 
ot ov I av 



(A.5) 



In spherical coordinates the integration of eq.(A.5) yields 
(iFridman fc Polyachenko (1984)1): 



/i = 



* d.'r^^+i(^„,)i^ 



dr dvr 



r 



dv± 



(A.6) 



where the integration is taken along the non-perturbed 
trajectory and we introduced the linear operator R defined 
by: 



R = cos a 



d sin a d 



89 sin t 



sin a cot 9 



d_ 
da 



(A.7) 



In eq.( [A.6| ) we used the fact that Jq{E, L) does not depend 
on a. Of course, <&i(?', 9, (p, t) does not depend on a either. 
Next, using: 

E=]-{vl+vl) + ^a{r), L = rv^, (A.8) 



as well as: 

d$i _ 9$i 

"dT ^ "dT 

along the particle trajectory, we can write eq.([A.6|) as 



or r 



(A.9) 



hit) ^ 



*'W-/'-'t 



(A.IO) 
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Thus, once we are given a perturbed gravitational poten- 
tial $i(r,t) the eq.(A.lC) yields the associated distribu- 
tion function /i(r,v,t). Next, to obtain the eigenmodes 
of the system vifc simply need to take into account the 
Poisson equation. Indeed, the potential $i(r, t) must obey 
the consistency requirement: 



This gives the action of the operator R on the functions 
T!^j „ through eq.(A.15) and eq.(A.16). In particular, we 
obtain: 



rt: 



' - ^^^^ (tLs + t^^-^) 



m,0 



2i 



(A.18) 



A$i = inGpi 



(A.ll) 



where the perturbed density pi is obtained from the distri- 
bution /i(r, V, i) derived in eq.( A.10| ). Thus, substituting 
eq.(A.lO) into eq.(A.ll) yields a linear eigenvalue problem 
for the unknown function $i(r,t). This procedure gives 
the eigenmodes of the system we look for. 

The radi al-orb it instability will be produced by the last 
term in eq.( A.10| ). Indeed, when the typical angular mo- 
mentum fi of the particles becomes very small the deriva- 
tive dfo/dL gives rise to a factor 1/fj, which diverges in the 
radial limit. This is similar to the usual Jeans instability: 
in the limit /z — > the velocity dispersion in the transverse 
dir ection s vanishes. Now, we look for the eigenmodes of 
eq. ( A. 10 ) . Because of spherical symmetry the angular part 
of the perturbation can be decomposed over the spher- 
ical harmonics Y™{9^<j)). However, in order to simplify 
the analysis we use the three- index functions T^ „((^, 0, a) 
which provide a representation of the rotation of Euler 



Thus, eq.( A.lO ) reads: 






VW+^dfo 



2i 



dL 



f di'e"*x(T;ki+rm,-i)- (A.19) 

-/ — oo 



We look for very large growth rates {uj — > +oo) hence we 
expand eq.( A.19 ) over powers of l/oj. To do so, we use the 
relation (obtained from successive integrations by parts): 



dt' e'^*' 



fit') - — 



f{t)- 



f{t) , /"(t) 



(A.20) 



for arbitrary functions /(t) such that the integral con- 
verges. Thus, up to second-order over 1/w we obtain: 

fm-,^tdfo[l d 1 52 \ 



angles (0,6', a), see Fridman & Polyachenko (1984) and 
Vilenkin (1968). Thus, we look for a perturbed gravita- 
tional potential of the form: 



^^{r,t)^e'^'x{r)T^A^,e,a). 



(A.12) 



Note that T.^^{ip,0,a) does not depend on a for n = 
0. Moreover, the functions T' 



(A.21) 

The perturbed density pi (r) is related to the distribution 
function /i by: 

27r 



m,0 



the usual spherical harmonics since we have yj'"(0, 0) ex 
e'^'^'^'^T}^Q{(j>,6,a). In particular, they also form a com- 
plete orthogonal system of functions on the sphere. The 
reason why we use the functions T^ „ is that they provide 
a very convenient basis to write the action of the operator 
R. Indeed, let us define the operators H+ and _ff_ by: 



are closely related to PiW == / div / -^rj da/i(r,v). (A.22) 



Besides, the functions TL 



3, a) are of the form: 



y^^„(0,0,a) = e-'('"^+"")/^t.,„( 



cost 



(A.23) 



,„ „ are closel y related to the 



H+ = e-'" 
and: 



d 



1 d 



0— 
do sin dcj) ' da 



.d_ J_d__ _d_ 

do sin 86 da 



(A.13) 



(A.14) 



where the functions P/, 

usual associated Legendre polynomials ([Vilenkin (1968)[ ) 
Therefore, only the terms T^ q (i.e. n = 0) contribute to 
the density pi after integration over the angle a. This is 
the reason why we needed to go up to second-order over 
l/u in eq.(A.21) since the first-order term which is pro- 



The n, the action of t hese operators H+ and H- is simply 
(see [Vilenkin (1968)1): 



portional to dfo/dL does not contain any factor T^ q. 

Now, we must examine under which conditions the 
density pi is indeed dominated by the last term written 



ff+T, 



and: 



VG-")a + "+l)Tm,n+l 



H.T^^^ = -V('+^)(^-"+i)7;U-i- 



(A.15) 



(A.16) 



Ne xt, th e comparison of eq.(A/7) with eq.( A.13 ) and 
eq.( A.14 ) yields the relation: 



in the expression (A.21). To do so, we must first recall the 
properties of the particle trajectories in a static spherical 
gravitational potential ^o{r). As is well-known, the orbit 
of a given particle is actually restricted to a plane which 
contains the origin (i.e. the center of the halo at r = 0) and 
which is orthogonal to the angular momentum L. Besides, 
the motion within this orbital plane can be described as a 
radial oscillation of frequency fi^ coupled with an angular 
oscillation of frequency Qe , with: 



R^-AH^ + H^ 

2,1 



(A.17) 



TT 



dr 



^2{E-^o{r))-L^/r^ 



(A.24) 
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and: 



Ldr 



r-^^2{E-^o{r))-L^/r^ 



(A.25) 



Here r,nin and r^ax are the minimum and maximum radii 
on the orbit. For extended halos we usually have 1/2 < 
Qg/flr < 1 since for a point mass (i.e. a Keplerian grav- 
itational potential) we get fig — Qr while for a constant 
density (i.e. the harmonic oscillator) we have fig = r2m/2 
(see Binney & Tremaine (1987)). Besides, for almost ra- 
dial orbits (i.e. L ^ 0) we usually have the resonance 
fig — flr/2 (the trajectory becomes almost symmetric 
with respect to the origin). For instance, we can check that 
this is the case for the density profile pQ{r) ex r~^ which 
would arise from eg. (|80|) . This resonance leads to the 
usual radial-orbit instability investigated for instance in 
Palmer fc Papaloizou (1987) and Polyachenko (1992) (see 
also references therein). However, the instability we study 
here is rather different as it does not rely on this resonance. 
Indeed, these works considered the secular instability of 
nearly radial orbits within a halo where most particles 
have a significant angular momentum. In other words, the 
mass of the "active" particles (i.e. those with almost ra- 
dial orbits which drive the instability) is small compared 
with the mass of the "passive" halo which determines the 
potential $0 (i.e. the particles with significant angular mo- 
mentum which are stable). In such a case, the radial-orbit 
instability only involves a small fraction of the matter and 
it exhibits a slow growth rate w, that is w ^ fio where 
J7o is the typical frequency of the system (i.e. f2o ^ \JQpc 



while the typical dynamical time is to 



where 



Pc is the average density of the halo). Indeed, the insta- 
bility is driven by the resonance fig = flr/2 so that the 
perturbation slowly increases through small cumulative ef- 
fects over many orbital periods. By contrast, the system 
we investigate here is a halo where all particles are "ac- 
tive" and follow nearly radial orbits. Of course, this leads 
to a much more violent instability. In particular, we shall 
see below that we now obtain a very large growth rate 
w 3> r^Q. This implies that the perturbation exhibits a 
strong growth in less than an orbital period. Therefore, 
the instability is not due to a resonance. In this sense, it is 
somewhat simpler and closer to the usual Jeans instability: 
it merely expresses the fact that in the absence of angular 
momentum there is no transverse velocity dispersion to 
stabilize the system against non-spherical perturbations. 
From the properties of the nearly radial orbits we can 
obtain the magnitude of the various terms in eq.( A.2l|) . 
Firstly, each derivative d/dt (along the particle trajectory) 
yields a factor flo, where we note flo the typical frequency 
of the system: 



flo 



(A.26) 



Indeed, we noticed above that for radial orbits fig — flr/2 
so that both frequencies are of the same order and c lose 
to Qq. Therefore, in order to use the expansion ( A.20| ) we 
must have flo/uj <ti 1. Secondly, the second-order term 



in dfo/dL dominates over the first-order term in dfo/dE 
if Eq ^ ujp. Here we note Eq the typical energy of the 
particles which is of order Eo ^ (Rflo)^, where R is the 
radius of the halo. These two constraints yield: 



flo '^ u} '^ — flo with in 

This also implies: 
p. < Lo- 



R'^flo. 



(A.27) 



(A.28) 



Note that Lo is the typical angular momentum of generic 
orbits with v± ^ Vr. Thus, we see at once from eq.(A.28) 
that the radial-orbit instability we shall obtain below will 
stop when the particles exhibit significant transverse ve- 
locities. Therefore, the main effect of this process is to 
isotropize the velocity dispersion. In this respect, it is sim- 
ilar to violent relaxation: starting from a distribution func- 
tion /o which is very far from thermodynamical equilib- 
rium (since the transverse velocity dispersion a± is zero) 
one obtains in a very short time (an orbital period !) a 
relaxed distribution where a± is of the order of the radial 
velocity dispersion ar- 

In the limit of nearly radial orbits {p <C Lq) we can 



look for eigenmodes which obey the constraint ( A.27 ). 
Then, the density pi is dominated by the last term in 
eq.( A.2l| ) which yields: 



Pi 



x/lil + 1) e"* f dw^dL^da dfo d 



2i 



nl 



2r2 



dL^dV^"'-^ 



T 



m,— 1 7 ■ 

(A.29) 

Thus, we now need the time derivative 9T/^ n/^^- From 
the analysis of the trajectory in the orbital plane one can 
show that: 
d 



(A.30) 



where the operator R was defined in cq.( A.7) and fig is the 
angular frequency of the orbit. Besides, using cq.(A.15) 
through eq.(A.17) we have: 



JdaRiT^^,+T^^_,) 



vw+^ 



da T, 



(A.31) 



where we only need to keep the term n = since we 
integrate over a. This yields: 



p, = l{l + 1)^ — X Tl^o 



e 



dvrdL'^ dfo 



2r2 dL 



-^ Qg. (A.32) 



As stated above, we can check in eq.( A.32) that the spher- 
ical harmonics separate. Indeed, starting from a potential 
$1 of the form (A.12|) we obtain a density pi which is 



proportional to the same harmonic r,'„ q. Next, substitut- 



ing the distribution function /o written in eq.(A.l) with 
eq.(A.3) we obtain: 



pi = -27tNo 



1(1 + 1) 



r^p LO^ ^ ™'° 
*o(fl) dS 



$o(r) y/2{E~^o{r)) 



fig g{E) 



(A.33) 
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where Nq is a positive number of order unity given by: vifhere we introduced the linear operator: 



No 



dx hix), No > 0. 



(A.34) 



In eq.(A.33) we used the fact that h{L) is strongly peaked 
for small values of L (of order of /x which obeys eq.( A.28)) 
so that f2e and \vr\ in eg. ( A. 84 ) are taken as functions of 
E with L = (i.e. along radial orbits). Finally, we obtain 
the eigenmodes through the Poisson equation ( A.ll). This 
yields: 



d 
dr 



r^^)-lil + l)x = -&7r^GNo 



dr J 

*o(-R.) 



/(/ + !) 



dE 



$o(r) ^2{E~<^o{r)) 



ne g{E) 



(A.35) 



which is a linear eigenvalue problem for the function of 
one variable x(?'). Thus, let us define the dimensionless 
linear operators Li and L2 by: 



{Ux){r 



and: 



dr V dr ' 



l)x(r) 



{L2.x){r) = -l{l + l)ar)x(r) 



(A.36) 



(A.37) 



where we introduced the dimensionless function Q{r) given 
by: 



C(0 



for r < R and C,{r) = 
equation ( A.35 ) reads: 



dE 



Vie g{E) (A.38) 



V2(£;~$o(r)) 
= for r > R. Then, the linear 



Li-X 



no 



L2-X 



(A.39) 



Next, it is easy to see from eq.(A.36) and eq.(A.37) that 
for I > 1 both linear operators Li and L2 are self-adjoint 
and negative definite over all functions x(j') which are not 
identically zero over < r < R, where we defined the 
scalar product: 



(Xi|i»IX2) = / dr xi{r){A-X2){r) 



(A.40) 



for any real functions xi and X2- Indeed, the function C(r) 
defined in eq.( A.38|) obeys C{r) > for r < i? and ({r) — 
for r > R. Therefore, if the fu nction x(r) is a solution of 
the eigenvalue problem ( A.39| ) we see that its associated 
growth rate ui obeys w^ > 0, by taking the scalar product 
(x|ii|x) in eq.( A.39 ). Indeed, if x(r) is zero over < 
r < R, then we have Li.x — over [0, +00 [ which implies 
X = (using the boundary conditions x(^) ~^ foi' r ^ 00 
and x(^ = 0) finite). Thus, we get a real growth rate w > 
which implies that the system is unstable for non-spherical 
perturbations (i.e. I > 1). Finally, the eigenmodes may be 
obtained from the more usual eigenvalue problem: 



L.x 



n^ 



(A.41) 



. ^ 1 . 



1 



C(r) 



1- 



1 



/(/ + l)dr 



dr 



(A.42) 



Thus, one first solves eq.( A.41 ) over < r < i?, which 
yields a continuous spectrum of functions x(^)i and next 
the match at R with the standard decaying solution of 
Li-x = defined over r > R gives rise to a discrete spec- 
trum. Let us note A^ the discrete eigenvalues of the op- 
erator L obtained in this way. As shown above, they are 
positive dimensionless numbers which give the associated 
growth rates through: 



V Afc/i 



(A.43) 



Thus, we see that the growth rates are of order uj ^ 
flo\/Lo/fJ-. Therefore, they satisfy the condition ( A. 27 ) 
when /x <C io, which validates the calculation described 
above. Note that perturbation theory is singular for a sys- 
tem with exactly radial orbits since the growth rates cu 
diverge in the limit /x — > 0. This is why we introduced a 
small angular momentum /i. On the other hand, we stress 
that the analysis described above is quite general. Thus, 
it does not rely on the shape of the equilibrium density 
profile po nor on the distribution function /q. Indeed, it is 
only based on the conditions (A. 27) and (A.28). Therefore, 
the system remains unstable until the constraint (A.28) is 
violated. This means that the halo quickly reaches an equi- 
librium state where the radial and transverse velocities are 
of the same order. 



References 

Balian, R., Schaefler, R., 1989, A&A 220, 1 

Bernardeau, P., 1994, ApJ 427, 51 

Bertschinger, E., 1985, ApJS 58, 39 

Binney, J., Tremaine, S., 1987, Galactic dynamics, Princeton 

University Press 
Cole, S., Lacey, C, 1996, MNRAS 281, 716 
Colombi, S., Bouchet, F.R., Hernquist, L., 1996, ApJ 465, 14 
Fillmore, J.A., Goldreich, P., 1984, ApJ 281, 1 
Fridman, A.M., Polyachenko, V.L., 1984, Physics of gravitating 

systems, Springer- Verlag 
Governato, F., Babul, A., Quinn, T., Baugh, C.M., Katz, N., 

Lake, G., 1999, MNRAS 307, 949 
Gradshteyn, I.S., Ryzhik, I.M., 1965, Table of integrals, series 

and products, fourth edition. Academic Press 
Hamilton, A.J.S., Kumar, P., Lu, E., Matthews, A., 1991, ApJ 

374, LI 
Jain, B., Mo, H.J., White, S.D.M., 1995, MNRAS 276, L25 
Kandrup, H.E., Sygnet, J. P., 1985, ApJ 298, 27 
Kandrup, H.E., Mahon, M.E., Smith, H., 1993, A&A 271, 440 
Palmer, P.L., Papaloizou, J., 1987, MNRAS 224, 1043 
Peacock, J.A., Dodds, S.J., 1996, MNRAS 280, L19 
Peebles, P.J.E., 1980, The large scale structure of the universe, 

Princeton University Press 
Polyachenko, V.L., 1992, Sov. Phys. JETP 74 (5), 755 
Press, W.H., Schechter, P., 1974, ApJ 187, 425 
Scoccimarro, R., Frieman, J.A., 1996, ApJ 473, 620 



p. Valageas: Dynamics of gravitational clustering IV. The probability distribution of rare events. 



27 



Teyssier, R., Chieze, J.-P., Alimi, J.-M., 1997, ApJ 480, 36 
Tormen, G., Bouchet, F.R., White, S.D.M., 1997, MNRAS 286, 

865 
Valageas, P., Schaeffer, R., 1997, A&A 328, 435 
Valageas, P., 1998, A&A 337, 655 

Valageas, P., Lacey, C, Schaeffer, R., 2000, MNRAS 311, 234 
Valageas, P., 2001, A&A 379, 8, Paper I 



Valage as, P., 2001, accepted by A&A, Paper II | astro 
ph/|0107126t 



Valage as, P., 2 001, accepted by A&A, Paper III | astro 
ph/|0107196t 



Valage as, P., 2001, accepted by A&A, Paper V astro- 
ph/|0109408| 



Vilenkin, N.J., 1968, Special functions and the theory of 
group representations, American mathematical society. 
Providence, Rhode Island 



